1 /* 2 Code for time stepping with the Runge-Kutta method(single-rate RK and multi-rate RK) 3 4 Notes: 5 1) The general system is written as 6 Udot = F(t,U) for single-rate RK; 7 2) The general system is written as 8 Udot = F(t,U) for combined RHS multi-rate RK, 9 user should give the indexes for both slow and fast components; 10 3) The general system is written as 11 Usdot = Fs(t,Us,Uf) 12 Ufdot = Ff(t,Us,Uf) for partitioned RHS multi-rate RK, 13 user should partioned RHS by themselves and also provide the indexes for both slow and fast components. 14 */ 15 16 #include <petsc/private/tsimpl.h> /*I "petscts.h" I*/ 17 #include <petscdm.h> 18 19 static TSRKType TSRKDefault = TSRK3BS; 20 static TSRKType TSMRKDefault = TSRK2A; 21 static PetscBool TSRKRegisterAllCalled; 22 static PetscBool TSRKPackageInitialized; 23 24 typedef struct _RKTableau *RKTableau; 25 struct _RKTableau { 26 char *name; 27 PetscInt order; /* Classical approximation order of the method i */ 28 PetscInt s; /* Number of stages */ 29 PetscInt p; /* Interpolation order */ 30 PetscBool FSAL; /* flag to indicate if tableau is FSAL */ 31 PetscReal *A,*b,*c; /* Tableau */ 32 PetscReal *bembed; /* Embedded formula of order one less (order-1) */ 33 PetscReal *binterp; /* Dense output formula */ 34 PetscReal ccfl; /* Placeholder for CFL coefficient relative to forward Euler */ 35 }; 36 typedef struct _RKTableauLink *RKTableauLink; 37 struct _RKTableauLink { 38 struct _RKTableau tab; 39 RKTableauLink next; 40 }; 41 static RKTableauLink RKTableauList; 42 43 typedef struct { 44 RKTableau tableau; 45 Vec *Y; /* States computed during the step */ 46 Vec *YdotRHS; /* Function evaluations for the non-stiff part and contains all components */ 47 Vec *YdotRHS_fast; /* Function evaluations for the non-stiff part and contains fast components */ 48 Vec *YdotRHS_slow; /* Function evaluations for the non-stiff part and contains slow components */ 49 Vec *VecDeltaLam; /* Increment of the adjoint sensitivity w.r.t IC at stage */ 50 Vec *VecDeltaMu; /* Increment of the adjoint sensitivity w.r.t P at stage */ 51 Vec VecCostIntegral0; /* backup for roll-backs due to events */ 52 PetscScalar *work; /* Scalar work */ 53 PetscInt slow; /* flag indicates call slow components solver (0) or fast components solver (1) */ 54 PetscReal stage_time; 55 TSStepStatus status; 56 PetscReal ptime; 57 PetscReal time_step; 58 PetscInt dtratio; /* ratio between slow time step size and fast step size */ 59 } TS_RK; 60 61 62 /*MC 63 TSRK1FE - First order forward Euler scheme. 64 65 This method has one stage. 66 67 Options database: 68 . -ts_rk_type 1fe 69 70 Level: advanced 71 72 .seealso: TSRK, TSRKType, TSRKSetType() 73 M*/ 74 /*MC 75 TSRK2A - Second order RK scheme. 76 77 This method has two stages. 78 79 Options database: 80 . -ts_rk_type 2a 81 82 Level: advanced 83 84 .seealso: TSRK, TSRKType, TSRKSetType() 85 M*/ 86 /*MC 87 TSRK3 - Third order RK scheme. 88 89 This method has three stages. 90 91 Options database: 92 . -ts_rk_type 3 93 94 Level: advanced 95 96 .seealso: TSRK, TSRKType, TSRKSetType() 97 M*/ 98 /*MC 99 TSRK3BS - Third order RK scheme of Bogacki-Shampine with 2nd order embedded method. 100 101 This method has four stages with the First Same As Last (FSAL) property. 102 103 Options database: 104 . -ts_rk_type 3bs 105 106 Level: advanced 107 108 References: https://doi.org/10.1016/0893-9659(89)90079-7 109 110 .seealso: TSRK, TSRKType, TSRKSetType() 111 M*/ 112 /*MC 113 TSRK4 - Fourth order RK scheme. 114 115 This is the classical Runge-Kutta method with four stages. 116 117 Options database: 118 . -ts_rk_type 4 119 120 Level: advanced 121 122 .seealso: TSRK, TSRKType, TSRKSetType() 123 M*/ 124 /*MC 125 TSRK5F - Fifth order Fehlberg RK scheme with a 4th order embedded method. 126 127 This method has six stages. 128 129 Options database: 130 . -ts_rk_type 5f 131 132 Level: advanced 133 134 .seealso: TSRK, TSRKType, TSRKSetType() 135 M*/ 136 /*MC 137 TSRK5DP - Fifth order Dormand-Prince RK scheme with the 4th order embedded method. 138 139 This method has seven stages with the First Same As Last (FSAL) property. 140 141 Options database: 142 . -ts_rk_type 5dp 143 144 Level: advanced 145 146 References: https://doi.org/10.1016/0771-050X(80)90013-3 147 148 .seealso: TSRK, TSRKType, TSRKSetType() 149 M*/ 150 /*MC 151 TSRK5BS - Fifth order Bogacki-Shampine RK scheme with 4th order embedded method. 152 153 This method has eight stages with the First Same As Last (FSAL) property. 154 155 Options database: 156 . -ts_rk_type 5bs 157 158 Level: advanced 159 160 References: https://doi.org/10.1016/0898-1221(96)00141-1 161 162 .seealso: TSRK, TSRKType, TSRKSetType() 163 M*/ 164 165 /*@C 166 TSRKRegisterAll - Registers all of the Runge-Kutta explicit methods in TSRK 167 168 Not Collective, but should be called by all processes which will need the schemes to be registered 169 170 Level: advanced 171 172 .keywords: TS, TSRK, register, all 173 174 .seealso: TSRKRegisterDestroy() 175 @*/ 176 PetscErrorCode TSRKRegisterAll(void) 177 { 178 PetscErrorCode ierr; 179 180 PetscFunctionBegin; 181 if (TSRKRegisterAllCalled) PetscFunctionReturn(0); 182 TSRKRegisterAllCalled = PETSC_TRUE; 183 184 #define RC PetscRealConstant 185 { 186 const PetscReal 187 A[1][1] = {{0}}, 188 b[1] = {RC(1.0)}; 189 ierr = TSRKRegister(TSRK1FE,1,1,&A[0][0],b,NULL,NULL,0,NULL);CHKERRQ(ierr); 190 } 191 { 192 const PetscReal 193 A[2][2] = {{0,0}, 194 {RC(1.0),0}}, 195 b[2] = {RC(0.5),RC(0.5)}, 196 bembed[2] = {RC(1.0),0}; 197 ierr = TSRKRegister(TSRK2A,2,2,&A[0][0],b,NULL,bembed,0,NULL);CHKERRQ(ierr); 198 } 199 { 200 const PetscReal 201 A[3][3] = {{0,0,0}, 202 {RC(2.0)/RC(3.0),0,0}, 203 {RC(-1.0)/RC(3.0),RC(1.0),0}}, 204 b[3] = {RC(0.25),RC(0.5),RC(0.25)}; 205 ierr = TSRKRegister(TSRK3,3,3,&A[0][0],b,NULL,NULL,0,NULL);CHKERRQ(ierr); 206 } 207 { 208 const PetscReal 209 A[4][4] = {{0,0,0,0}, 210 {RC(1.0)/RC(2.0),0,0,0}, 211 {0,RC(3.0)/RC(4.0),0,0}, 212 {RC(2.0)/RC(9.0),RC(1.0)/RC(3.0),RC(4.0)/RC(9.0),0}}, 213 b[4] = {RC(2.0)/RC(9.0),RC(1.0)/RC(3.0),RC(4.0)/RC(9.0),0}, 214 bembed[4] = {RC(7.0)/RC(24.0),RC(1.0)/RC(4.0),RC(1.0)/RC(3.0),RC(1.0)/RC(8.0)}; 215 ierr = TSRKRegister(TSRK3BS,3,4,&A[0][0],b,NULL,bembed,0,NULL);CHKERRQ(ierr); 216 } 217 { 218 const PetscReal 219 A[4][4] = {{0,0,0,0}, 220 {RC(0.5),0,0,0}, 221 {0,RC(0.5),0,0}, 222 {0,0,RC(1.0),0}}, 223 b[4] = {RC(1.0)/RC(6.0),RC(1.0)/RC(3.0),RC(1.0)/RC(3.0),RC(1.0)/RC(6.0)}; 224 ierr = TSRKRegister(TSRK4,4,4,&A[0][0],b,NULL,NULL,0,NULL);CHKERRQ(ierr); 225 } 226 { 227 const PetscReal 228 A[6][6] = {{0,0,0,0,0,0}, 229 {RC(0.25),0,0,0,0,0}, 230 {RC(3.0)/RC(32.0),RC(9.0)/RC(32.0),0,0,0,0}, 231 {RC(1932.0)/RC(2197.0),RC(-7200.0)/RC(2197.0),RC(7296.0)/RC(2197.0),0,0,0}, 232 {RC(439.0)/RC(216.0),RC(-8.0),RC(3680.0)/RC(513.0),RC(-845.0)/RC(4104.0),0,0}, 233 {RC(-8.0)/RC(27.0),RC(2.0),RC(-3544.0)/RC(2565.0),RC(1859.0)/RC(4104.0),RC(-11.0)/RC(40.0),0}}, 234 b[6] = {RC(16.0)/RC(135.0),0,RC(6656.0)/RC(12825.0),RC(28561.0)/RC(56430.0),RC(-9.0)/RC(50.0),RC(2.0)/RC(55.0)}, 235 bembed[6] = {RC(25.0)/RC(216.0),0,RC(1408.0)/RC(2565.0),RC(2197.0)/RC(4104.0),RC(-1.0)/RC(5.0),0}; 236 ierr = TSRKRegister(TSRK5F,5,6,&A[0][0],b,NULL,bembed,0,NULL);CHKERRQ(ierr); 237 } 238 { 239 const PetscReal 240 A[7][7] = {{0,0,0,0,0,0,0}, 241 {RC(1.0)/RC(5.0),0,0,0,0,0,0}, 242 {RC(3.0)/RC(40.0),RC(9.0)/RC(40.0),0,0,0,0,0}, 243 {RC(44.0)/RC(45.0),RC(-56.0)/RC(15.0),RC(32.0)/RC(9.0),0,0,0,0}, 244 {RC(19372.0)/RC(6561.0),RC(-25360.0)/RC(2187.0),RC(64448.0)/RC(6561.0),RC(-212.0)/RC(729.0),0,0,0}, 245 {RC(9017.0)/RC(3168.0),RC(-355.0)/RC(33.0),RC(46732.0)/RC(5247.0),RC(49.0)/RC(176.0),RC(-5103.0)/RC(18656.0),0,0}, 246 {RC(35.0)/RC(384.0),0,RC(500.0)/RC(1113.0),RC(125.0)/RC(192.0),RC(-2187.0)/RC(6784.0),RC(11.0)/RC(84.0),0}}, 247 b[7] = {RC(35.0)/RC(384.0),0,RC(500.0)/RC(1113.0),RC(125.0)/RC(192.0),RC(-2187.0)/RC(6784.0),RC(11.0)/RC(84.0),0}, 248 bembed[7] = {RC(5179.0)/RC(57600.0),0,RC(7571.0)/RC(16695.0),RC(393.0)/RC(640.0),RC(-92097.0)/RC(339200.0),RC(187.0)/RC(2100.0),RC(1.0)/RC(40.0)}; 249 ierr = TSRKRegister(TSRK5DP,5,7,&A[0][0],b,NULL,bembed,0,NULL);CHKERRQ(ierr); 250 } 251 { 252 const PetscReal 253 A[8][8] = {{0,0,0,0,0,0,0,0}, 254 {RC(1.0)/RC(6.0),0,0,0,0,0,0,0}, 255 {RC(2.0)/RC(27.0),RC(4.0)/RC(27.0),0,0,0,0,0,0}, 256 {RC(183.0)/RC(1372.0),RC(-162.0)/RC(343.0),RC(1053.0)/RC(1372.0),0,0,0,0,0}, 257 {RC(68.0)/RC(297.0),RC(-4.0)/RC(11.0),RC(42.0)/RC(143.0),RC(1960.0)/RC(3861.0),0,0,0,0}, 258 {RC(597.0)/RC(22528.0),RC(81.0)/RC(352.0),RC(63099.0)/RC(585728.0),RC(58653.0)/RC(366080.0),RC(4617.0)/RC(20480.0),0,0,0}, 259 {RC(174197.0)/RC(959244.0),RC(-30942.0)/RC(79937.0),RC(8152137.0)/RC(19744439.0),RC(666106.0)/RC(1039181.0),RC(-29421.0)/RC(29068.0),RC(482048.0)/RC(414219.0),0,0}, 260 {RC(587.0)/RC(8064.0),0,RC(4440339.0)/RC(15491840.0),RC(24353.0)/RC(124800.0),RC(387.0)/RC(44800.0),RC(2152.0)/RC(5985.0),RC(7267.0)/RC(94080.0),0}}, 261 b[8] = {RC(587.0)/RC(8064.0),0,RC(4440339.0)/RC(15491840.0),RC(24353.0)/RC(124800.0),RC(387.0)/RC(44800.0),RC(2152.0)/RC(5985.0),RC(7267.0)/RC(94080.0),0}, 262 bembed[8] = {RC(2479.0)/RC(34992.0),0,RC(123.0)/RC(416.0),RC(612941.0)/RC(3411720.0),RC(43.0)/RC(1440.0),RC(2272.0)/RC(6561.0),RC(79937.0)/RC(1113912.0),RC(3293.0)/RC(556956.0)}; 263 ierr = TSRKRegister(TSRK5BS,5,8,&A[0][0],b,NULL,bembed,0,NULL);CHKERRQ(ierr); 264 } 265 #undef RC 266 PetscFunctionReturn(0); 267 } 268 269 /*@C 270 TSRKRegisterDestroy - Frees the list of schemes that were registered by TSRKRegister(). 271 272 Not Collective 273 274 Level: advanced 275 276 .keywords: TSRK, register, destroy 277 .seealso: TSRKRegister(), TSRKRegisterAll() 278 @*/ 279 PetscErrorCode TSRKRegisterDestroy(void) 280 { 281 PetscErrorCode ierr; 282 RKTableauLink link; 283 284 PetscFunctionBegin; 285 while ((link = RKTableauList)) { 286 RKTableau t = &link->tab; 287 RKTableauList = link->next; 288 ierr = PetscFree3(t->A,t->b,t->c); CHKERRQ(ierr); 289 ierr = PetscFree (t->bembed); CHKERRQ(ierr); 290 ierr = PetscFree (t->binterp); CHKERRQ(ierr); 291 ierr = PetscFree (t->name); CHKERRQ(ierr); 292 ierr = PetscFree (link); CHKERRQ(ierr); 293 } 294 TSRKRegisterAllCalled = PETSC_FALSE; 295 PetscFunctionReturn(0); 296 } 297 298 /*@C 299 TSRKInitializePackage - This function initializes everything in the TSRK package. It is called 300 from TSInitializePackage(). 301 302 Level: developer 303 304 .keywords: TS, TSRK, initialize, package 305 .seealso: PetscInitialize() 306 @*/ 307 PetscErrorCode TSRKInitializePackage(void) 308 { 309 PetscErrorCode ierr; 310 311 PetscFunctionBegin; 312 if (TSRKPackageInitialized) PetscFunctionReturn(0); 313 TSRKPackageInitialized = PETSC_TRUE; 314 ierr = TSRKRegisterAll();CHKERRQ(ierr); 315 ierr = PetscRegisterFinalize(TSRKFinalizePackage);CHKERRQ(ierr); 316 PetscFunctionReturn(0); 317 } 318 319 /*@C 320 TSRKFinalizePackage - This function destroys everything in the TSRK package. It is 321 called from PetscFinalize(). 322 323 Level: developer 324 325 .keywords: Petsc, destroy, package 326 .seealso: PetscFinalize() 327 @*/ 328 PetscErrorCode TSRKFinalizePackage(void) 329 { 330 PetscErrorCode ierr; 331 332 PetscFunctionBegin; 333 TSRKPackageInitialized = PETSC_FALSE; 334 ierr = TSRKRegisterDestroy();CHKERRQ(ierr); 335 PetscFunctionReturn(0); 336 } 337 338 /*@C 339 TSRKRegister - register an RK scheme by providing the entries in the Butcher tableau and optionally embedded approximations and interpolation 340 341 Not Collective, but the same schemes should be registered on all processes on which they will be used 342 343 Input Parameters: 344 + name - identifier for method 345 . order - approximation order of method 346 . s - number of stages, this is the dimension of the matrices below 347 . A - stage coefficients (dimension s*s, row-major) 348 . b - step completion table (dimension s; NULL to use last row of A) 349 . c - abscissa (dimension s; NULL to use row sums of A) 350 . bembed - completion table for embedded method (dimension s; NULL if not available) 351 . p - Order of the interpolation scheme, equal to the number of columns of binterp 352 - binterp - Coefficients of the interpolation formula (dimension s*p; NULL to reuse b with p=1) 353 354 Notes: 355 Several RK methods are provided, this function is only needed to create new methods. 356 357 Level: advanced 358 359 .keywords: TS, register 360 361 .seealso: TSRK 362 @*/ 363 PetscErrorCode TSRKRegister(TSRKType name,PetscInt order,PetscInt s, 364 const PetscReal A[],const PetscReal b[],const PetscReal c[], 365 const PetscReal bembed[],PetscInt p,const PetscReal binterp[]) 366 { 367 PetscErrorCode ierr; 368 RKTableauLink link; 369 RKTableau t; 370 PetscInt i,j; 371 372 PetscFunctionBegin; 373 PetscValidCharPointer(name,1); 374 PetscValidRealPointer(A,4); 375 if (b) PetscValidRealPointer(b,5); 376 if (c) PetscValidRealPointer(c,6); 377 if (bembed) PetscValidRealPointer(bembed,7); 378 if (binterp || p > 1) PetscValidRealPointer(binterp,9); 379 380 ierr = TSRKInitializePackage();CHKERRQ(ierr); 381 ierr = PetscNew(&link);CHKERRQ(ierr); 382 t = &link->tab; 383 384 ierr = PetscStrallocpy(name,&t->name);CHKERRQ(ierr); 385 t->order = order; 386 t->s = s; 387 ierr = PetscMalloc3(s*s,&t->A,s,&t->b,s,&t->c);CHKERRQ(ierr); 388 ierr = PetscMemcpy(t->A,A,s*s*sizeof(A[0]));CHKERRQ(ierr); 389 if (b) { ierr = PetscMemcpy(t->b,b,s*sizeof(b[0]));CHKERRQ(ierr); } 390 else for (i=0; i<s; i++) t->b[i] = A[(s-1)*s+i]; 391 if (c) { ierr = PetscMemcpy(t->c,c,s*sizeof(c[0]));CHKERRQ(ierr); } 392 else for (i=0; i<s; i++) for (j=0,t->c[i]=0; j<s; j++) t->c[i] += A[i*s+j]; 393 t->FSAL = PETSC_TRUE; 394 for (i=0; i<s; i++) if (t->A[(s-1)*s+i] != t->b[i]) t->FSAL = PETSC_FALSE; 395 396 if (bembed) { 397 ierr = PetscMalloc1(s,&t->bembed);CHKERRQ(ierr); 398 ierr = PetscMemcpy(t->bembed,bembed,s*sizeof(bembed[0]));CHKERRQ(ierr); 399 } 400 401 if (!binterp) { p = 1; binterp = t->b; } 402 t->p = p; 403 ierr = PetscMalloc1(s*p,&t->binterp);CHKERRQ(ierr); 404 ierr = PetscMemcpy(t->binterp,binterp,s*p*sizeof(binterp[0]));CHKERRQ(ierr); 405 406 link->next = RKTableauList; 407 RKTableauList = link; 408 PetscFunctionReturn(0); 409 } 410 411 /* 412 This is for single-step RK method 413 The step completion formula is 414 415 x1 = x0 + h b^T YdotRHS 416 417 This function can be called before or after ts->vec_sol has been updated. 418 Suppose we have a completion formula (b) and an embedded formula (be) of different order. 419 We can write 420 421 x1e = x0 + h be^T YdotRHS 422 = x1 - h b^T YdotRHS + h be^T YdotRHS 423 = x1 + h (be - b)^T YdotRHS 424 425 so we can evaluate the method with different order even after the step has been optimistically completed. 426 */ 427 static PetscErrorCode TSEvaluateStep_RK(TS ts,PetscInt order,Vec X,PetscBool *done) 428 { 429 TS_RK *rk = (TS_RK*)ts->data; 430 RKTableau tab = rk->tableau; 431 PetscScalar *w = rk->work; 432 PetscReal h; 433 PetscInt s = tab->s,j; 434 PetscErrorCode ierr; 435 436 PetscFunctionBegin; 437 switch (rk->status) { 438 case TS_STEP_INCOMPLETE: 439 case TS_STEP_PENDING: 440 h = ts->time_step; break; 441 case TS_STEP_COMPLETE: 442 h = ts->ptime - ts->ptime_prev; break; 443 default: SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_PLIB,"Invalid TSStepStatus"); 444 } 445 if (order == tab->order) { 446 if (rk->status == TS_STEP_INCOMPLETE) { 447 ierr = VecCopy(ts->vec_sol,X);CHKERRQ(ierr); 448 for (j=0; j<s; j++) w[j] = h*tab->b[j]; 449 ierr = VecMAXPY(X,s,w,rk->YdotRHS);CHKERRQ(ierr); 450 } else {ierr = VecCopy(ts->vec_sol,X);CHKERRQ(ierr);} 451 PetscFunctionReturn(0); 452 } else if (order == tab->order-1) { 453 if (!tab->bembed) goto unavailable; 454 if (rk->status == TS_STEP_INCOMPLETE) { /*Complete with the embedded method (be)*/ 455 ierr = VecCopy(ts->vec_sol,X);CHKERRQ(ierr); 456 for (j=0; j<s; j++) w[j] = h*tab->bembed[j]; 457 ierr = VecMAXPY(X,s,w,rk->YdotRHS);CHKERRQ(ierr); 458 } else { /*Rollback and re-complete using (be-b) */ 459 ierr = VecCopy(ts->vec_sol,X);CHKERRQ(ierr); 460 for (j=0; j<s; j++) w[j] = h*(tab->bembed[j] - tab->b[j]); 461 ierr = VecMAXPY(X,s,w,rk->YdotRHS);CHKERRQ(ierr); 462 } 463 if (done) *done = PETSC_TRUE; 464 PetscFunctionReturn(0); 465 } 466 unavailable: 467 if (done) *done = PETSC_FALSE; 468 else SETERRQ3(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"RK '%s' of order %D cannot evaluate step at order %D. Consider using -ts_adapt_type none or a different method that has an embedded estimate.",tab->name,tab->order,order); 469 PetscFunctionReturn(0); 470 } 471 472 static PetscErrorCode TSForwardCostIntegral_RK(TS ts) 473 { 474 TS_RK *rk = (TS_RK*)ts->data; 475 RKTableau tab = rk->tableau; 476 const PetscInt s = tab->s; 477 const PetscReal *b = tab->b,*c = tab->c; 478 Vec *Y = rk->Y; 479 PetscInt i; 480 PetscErrorCode ierr; 481 482 PetscFunctionBegin; 483 /* backup cost integral */ 484 ierr = VecCopy(ts->vec_costintegral,rk->VecCostIntegral0);CHKERRQ(ierr); 485 for (i=s-1; i>=0; i--) { 486 /* Evolve ts->vec_costintegral to compute integrals */ 487 ierr = TSComputeCostIntegrand(ts,rk->ptime+rk->time_step*c[i],Y[i],ts->vec_costintegrand);CHKERRQ(ierr); 488 ierr = VecAXPY(ts->vec_costintegral,rk->time_step*b[i],ts->vec_costintegrand);CHKERRQ(ierr); 489 } 490 PetscFunctionReturn(0); 491 } 492 493 static PetscErrorCode TSAdjointCostIntegral_RK(TS ts) 494 { 495 TS_RK *rk = (TS_RK*)ts->data; 496 RKTableau tab = rk->tableau; 497 const PetscInt s = tab->s; 498 const PetscReal *b = tab->b,*c = tab->c; 499 Vec *Y = rk->Y; 500 PetscInt i; 501 PetscErrorCode ierr; 502 503 PetscFunctionBegin; 504 for (i=s-1; i>=0; i--) { 505 /* Evolve ts->vec_costintegral to compute integrals */ 506 ierr = TSComputeCostIntegrand(ts,ts->ptime+ts->time_step*(1.0-c[i]),Y[i],ts->vec_costintegrand);CHKERRQ(ierr); 507 ierr = VecAXPY(ts->vec_costintegral,-ts->time_step*b[i],ts->vec_costintegrand);CHKERRQ(ierr); 508 } 509 PetscFunctionReturn(0); 510 } 511 512 static PetscErrorCode TSRollBack_RK(TS ts) 513 { 514 TS_RK *rk = (TS_RK*)ts->data; 515 RKTableau tab = rk->tableau; 516 const PetscInt s = tab->s; 517 const PetscReal *b = tab->b; 518 PetscScalar *w = rk->work; 519 Vec *YdotRHS = rk->YdotRHS; 520 PetscInt j; 521 PetscReal h; 522 PetscErrorCode ierr; 523 524 PetscFunctionBegin; 525 switch (rk->status) { 526 case TS_STEP_INCOMPLETE: 527 case TS_STEP_PENDING: 528 h = ts->time_step; break; 529 case TS_STEP_COMPLETE: 530 h = ts->ptime - ts->ptime_prev; break; 531 default: SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_PLIB,"Invalid TSStepStatus"); 532 } 533 for (j=0; j<s; j++) w[j] = -h*b[j]; 534 ierr = VecMAXPY(ts->vec_sol,s,w,YdotRHS);CHKERRQ(ierr); 535 PetscFunctionReturn(0); 536 } 537 538 static PetscErrorCode TSStep_RK(TS ts) 539 { 540 TS_RK *rk = (TS_RK*)ts->data; 541 RKTableau tab = rk->tableau; 542 const PetscInt s = tab->s; 543 const PetscReal *A = tab->A,*c = tab->c; 544 PetscScalar *w = rk->work; 545 Vec *Y = rk->Y,*YdotRHS = rk->YdotRHS; 546 PetscBool FSAL = tab->FSAL; 547 TSAdapt adapt; 548 PetscInt i,j; 549 PetscInt rejections = 0; 550 PetscBool stageok,accept = PETSC_TRUE; 551 PetscReal next_time_step = ts->time_step; 552 PetscErrorCode ierr; 553 554 PetscFunctionBegin; 555 if (ts->steprollback || ts->steprestart) FSAL = PETSC_FALSE; 556 if (FSAL) { ierr = VecCopy(YdotRHS[s-1],YdotRHS[0]);CHKERRQ(ierr); } 557 558 rk->status = TS_STEP_INCOMPLETE; 559 while (!ts->reason && rk->status != TS_STEP_COMPLETE) { 560 PetscReal t = ts->ptime; 561 PetscReal h = ts->time_step; 562 for (i=0; i<s; i++) { 563 rk->stage_time = t + h*c[i]; 564 ierr = TSPreStage(ts,rk->stage_time); CHKERRQ(ierr); 565 ierr = VecCopy(ts->vec_sol,Y[i]);CHKERRQ(ierr); 566 for (j=0; j<i; j++) w[j] = h*A[i*s+j]; 567 ierr = VecMAXPY(Y[i],i,w,YdotRHS);CHKERRQ(ierr); 568 ierr = TSPostStage(ts,rk->stage_time,i,Y); CHKERRQ(ierr); 569 ierr = TSGetAdapt(ts,&adapt);CHKERRQ(ierr); 570 ierr = TSAdaptCheckStage(adapt,ts,rk->stage_time,Y[i],&stageok);CHKERRQ(ierr); 571 if (!stageok) goto reject_step; 572 if (FSAL && !i) continue; 573 ierr = TSComputeRHSFunction(ts,t+h*c[i],Y[i],YdotRHS[i]);CHKERRQ(ierr); 574 } 575 576 rk->status = TS_STEP_INCOMPLETE; 577 ierr = TSEvaluateStep(ts,tab->order,ts->vec_sol,NULL);CHKERRQ(ierr); 578 rk->status = TS_STEP_PENDING; 579 ierr = TSGetAdapt(ts,&adapt);CHKERRQ(ierr); 580 ierr = TSAdaptCandidatesClear(adapt);CHKERRQ(ierr); 581 ierr = TSAdaptCandidateAdd(adapt,tab->name,tab->order,1,tab->ccfl,(PetscReal)tab->s,PETSC_TRUE);CHKERRQ(ierr); 582 ierr = TSAdaptChoose(adapt,ts,ts->time_step,NULL,&next_time_step,&accept);CHKERRQ(ierr); 583 rk->status = accept ? TS_STEP_COMPLETE : TS_STEP_INCOMPLETE; 584 if (!accept) { /*Roll back the current step */ 585 ierr = TSRollBack_RK(ts);CHKERRQ(ierr); 586 ts->time_step = next_time_step; 587 goto reject_step; 588 } 589 590 if (ts->costintegralfwd) { /*Save the info for the later use in cost integral evaluation*/ 591 rk->ptime = ts->ptime; 592 rk->time_step = ts->time_step; 593 } 594 595 ts->ptime += ts->time_step; 596 ts->time_step = next_time_step; 597 break; 598 599 reject_step: 600 ts->reject++; accept = PETSC_FALSE; 601 if (!ts->reason && ++rejections > ts->max_reject && ts->max_reject >= 0) { 602 ts->reason = TS_DIVERGED_STEP_REJECTED; 603 ierr = PetscInfo2(ts,"Step=%D, step rejections %D greater than current TS allowed, stopping solve\n",ts->steps,rejections);CHKERRQ(ierr); 604 } 605 } 606 PetscFunctionReturn(0); 607 } 608 609 static PetscErrorCode TSInterpolate_RK(TS ts,PetscReal itime,Vec X) 610 { 611 TS_RK *rk = (TS_RK*)ts->data; 612 PetscInt s = rk->tableau->s,p = rk->tableau->p,i,j; 613 PetscReal h; 614 PetscReal tt,t; 615 PetscScalar *b; 616 const PetscReal *B = rk->tableau->binterp; 617 PetscErrorCode ierr; 618 619 PetscFunctionBegin; 620 if (!B) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"TSRK %s does not have an interpolation formula",rk->tableau->name); 621 622 switch (rk->status) { 623 case TS_STEP_INCOMPLETE: 624 case TS_STEP_PENDING: 625 h = ts->time_step; 626 t = (itime - ts->ptime)/h; 627 break; 628 case TS_STEP_COMPLETE: 629 h = ts->ptime - ts->ptime_prev; 630 t = (itime - ts->ptime)/h + 1; /* In the interval [0,1] */ 631 break; 632 default: SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_PLIB,"Invalid TSStepStatus"); 633 } 634 ierr = PetscMalloc1(s,&b);CHKERRQ(ierr); 635 for (i=0; i<s; i++) b[i] = 0; 636 for (j=0,tt=t; j<p; j++,tt*=t) { 637 for (i=0; i<s; i++) { 638 b[i] += h * B[i*p+j] * tt; 639 } 640 } 641 ierr = VecCopy(rk->Y[0],X);CHKERRQ(ierr); 642 ierr = VecMAXPY(X,s,b,rk->YdotRHS);CHKERRQ(ierr); 643 ierr = PetscFree(b);CHKERRQ(ierr); 644 PetscFunctionReturn(0); 645 } 646 647 /* 648 This is interpolate formula for partitioned RHS multirate RK method 649 */ 650 651 static PetscErrorCode TSInterpolate_PARTITIONEDMRK(TS ts,PetscReal itime,Vec X) 652 { 653 TS_RK *rk = (TS_RK*)ts->data; 654 PetscInt s = rk->tableau->s,p = rk->tableau->p,i,j; 655 Vec Yslow; /* vector holds the slow components in Y[0] */ 656 PetscReal h; 657 PetscReal tt,t; 658 PetscScalar *b; 659 const PetscReal *B = rk->tableau->binterp; 660 PetscErrorCode ierr; 661 662 PetscFunctionBegin; 663 if (!B) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"TSRK %s does not have an interpolation formula",rk->tableau->name); 664 665 switch (rk->status) { 666 case TS_STEP_INCOMPLETE: 667 case TS_STEP_PENDING: 668 h = ts->time_step; 669 t = (itime - ts->ptime)/h; 670 break; 671 case TS_STEP_COMPLETE: 672 h = ts->ptime - ts->ptime_prev; 673 t = (itime - ts->ptime)/h + 1; /* In the interval [0,1] */ 674 break; 675 default: SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_PLIB,"Invalid TSStepStatus"); 676 } 677 ierr = PetscMalloc1(s,&b);CHKERRQ(ierr); 678 for (i=0; i<s; i++) b[i] = 0; 679 for (j=0,tt=t; j<p; j++,tt*=t) { 680 for (i=0; i<s; i++) { 681 b[i] += h * B[i*p+j] * tt; 682 } 683 } 684 ierr = VecGetSubVector(rk->Y[0],ts->iss,&Yslow);CHKERRQ(ierr); 685 ierr = VecCopy(Yslow,X);CHKERRQ(ierr); 686 ierr = VecMAXPY(X,s,b,rk->YdotRHS_slow);CHKERRQ(ierr); 687 ierr = VecRestoreSubVector(rk->Y[0],ts->iss,&Yslow);CHKERRQ(ierr); 688 ierr = PetscFree(b);CHKERRQ(ierr); 689 PetscFunctionReturn(0); 690 } 691 692 /* 693 This is for combined RHS multirate RK method 694 The step completion formula is 695 696 x1 = x0 + h b^T YdotRHS 697 698 */ 699 static PetscErrorCode TSEvaluateStep_RKMULTIRATE(TS ts,PetscInt order,Vec X,PetscBool *done) 700 { 701 TS_RK *rk = (TS_RK*)ts->data; 702 RKTableau tab = rk->tableau; 703 Vec Xtemp; /* temporary work vector for X */ 704 PetscScalar *w = rk->work; 705 PetscScalar *x_ptr,*xtemp_ptr; /* location to put the pointer to X and Xtemp respectively */ 706 PetscReal h = ts->time_step; 707 PetscInt s = tab->s,j; 708 PetscInt len_slow,len_fast; /* the number of slow components and fast components respectively */ 709 const PetscInt *is_slow,*is_fast; /* the index for slow components and fast components respectively */ 710 PetscErrorCode ierr; 711 712 PetscFunctionBegin; 713 ierr = VecCopy(ts->vec_sol,X);CHKERRQ(ierr); 714 ierr = VecDuplicate(ts->vec_sol,&Xtemp);CHKERRQ(ierr); 715 ierr = VecCopy(ts->vec_sol,Xtemp);CHKERRQ(ierr); 716 if (!rk->slow){ 717 for (j=0; j<s; j++) w[j] = h*tab->b[j]; 718 /* update the value of slow components, and discard the updating value of fast components */ 719 ierr = VecMAXPY(Xtemp,s,w,rk->YdotRHS);CHKERRQ(ierr); 720 ierr = VecGetArray(X,&x_ptr);CHKERRQ(ierr); 721 ierr = VecGetArray(Xtemp,&xtemp_ptr);CHKERRQ(ierr); 722 ierr = ISGetSize(ts->iss,&len_slow);CHKERRQ(ierr); 723 ierr = ISGetIndices(ts->iss,&is_slow);CHKERRQ(ierr); 724 ierr = ISRestoreIndices(ts->iss,&is_slow);CHKERRQ(ierr); 725 for (j=0; j<len_slow; j++) x_ptr[is_slow[j]] = xtemp_ptr[is_slow[j]]; 726 ierr = VecRestoreArray(X,&x_ptr);CHKERRQ(ierr); 727 ierr = VecRestoreArray(Xtemp,&xtemp_ptr);CHKERRQ(ierr); 728 } else { 729 /* update the value of fast components, and discard the updating value of slow components */ 730 for (j=0; j<s; j++) w[j] = h/rk->dtratio*tab->b[j]; 731 ierr = VecMAXPY(Xtemp,s,w,rk->YdotRHS);CHKERRQ(ierr); 732 ierr = VecGetArray(X,&x_ptr);CHKERRQ(ierr); 733 ierr = VecGetArray(Xtemp,&xtemp_ptr);CHKERRQ(ierr); 734 ierr = ISGetSize(ts->isf,&len_fast);CHKERRQ(ierr); 735 ierr = ISGetIndices(ts->isf,&is_fast);CHKERRQ(ierr); 736 ierr = ISRestoreIndices(ts->isf,&is_fast);CHKERRQ(ierr); 737 for (j=0; j<len_fast; j++) x_ptr[is_fast[j]] = xtemp_ptr[is_fast[j]]; 738 ierr = VecRestoreArray(X,&x_ptr);CHKERRQ(ierr); 739 ierr = VecRestoreArray(Xtemp,&xtemp_ptr);CHKERRQ(ierr); 740 } 741 /* free temporary vector space */ 742 ierr = VecDestroy(&Xtemp);CHKERRQ(ierr); 743 PetscFunctionReturn(0); 744 } 745 746 static PetscErrorCode TSStep_RKMULTIRATE(TS ts) 747 { 748 TS_RK *rk = (TS_RK*)ts->data; 749 RKTableau tab = rk->tableau; 750 Vec *Y = rk->Y,*YdotRHS = rk->YdotRHS; 751 Vec stage_fast,stage_slow; /* vectors store the stage value of fast and slow components respectively */ 752 Vec presol_fast,newslo_fast; /* vectors store the previous and new time solution of fast components respectively */ 753 Vec YdotRHS_prev,prev_sol; /* vectors store the previous YdotRHS and solution respectively */ 754 const PetscInt s = tab->s,*is_slow,*is_fast; /* is_slow: index of slow components; is_fast: index of fast components */ 755 const PetscReal *A = tab->A,*c = tab->c; 756 PetscScalar *w = rk->work; 757 PetscScalar *y_ptr,*faststage_ptr,*slowstage_ptr; /* location to put the pointer to Y, stage_fast and stage_slow respectively */ 758 PetscScalar *ydotrhsprev_ptr,*ydotrhs_ptr; /* location to put the pointer to YdotRHS_prev and YdotRHS respectively */ 759 PetscInt i,j,k,len_slow,len_fast; /* len_slow: the number of slow comonents; len_fast: the number of fast components */ 760 PetscReal next_time_step = ts->time_step,t = ts->ptime,h = ts->time_step; 761 PetscErrorCode ierr; 762 763 PetscFunctionBegin; 764 rk->status = TS_STEP_INCOMPLETE; 765 ierr = VecDuplicate(ts->vec_sol,&stage_fast);CHKERRQ(ierr); 766 ierr = VecDuplicate(ts->vec_sol,&stage_slow);CHKERRQ(ierr); 767 ierr = VecDuplicate(ts->vec_sol,&YdotRHS_prev);CHKERRQ(ierr); 768 ierr = VecDuplicate(ts->vec_sol,&prev_sol);CHKERRQ(ierr); 769 ierr = ISGetSize(ts->iss,&len_slow);CHKERRQ(ierr); 770 ierr = ISGetSize(ts->isf,&len_fast);CHKERRQ(ierr); 771 ierr = ISGetIndices(ts->iss,&is_slow);CHKERRQ(ierr); 772 ierr = ISGetIndices(ts->isf,&is_fast);CHKERRQ(ierr); 773 ierr = ISRestoreIndices(ts->isf,&is_fast);CHKERRQ(ierr); 774 ierr = ISRestoreIndices(ts->iss,&is_slow);CHKERRQ(ierr); 775 for (i=0; i<s; i++) { 776 rk->stage_time = t + h*c[i]; 777 ierr = TSPreStage(ts,rk->stage_time);CHKERRQ(ierr); 778 ierr = VecCopy(ts->vec_sol,Y[i]);CHKERRQ(ierr); 779 for (j=0; j<i; j++) w[j] = h*A[i*s+j]; 780 ierr = VecMAXPY(Y[i],i,w,YdotRHS);CHKERRQ(ierr); 781 ierr = TSPostStage(ts,rk->stage_time,i,Y); CHKERRQ(ierr); 782 /* compute the stage RHS */ 783 ierr = TSComputeRHSFunction(ts,t+h*c[i],Y[i],YdotRHS[i]);CHKERRQ(ierr); 784 } 785 ierr = VecCopy(ts->vec_sol,prev_sol);CHKERRQ(ierr); 786 rk->slow = 0; 787 ierr = TSEvaluateStep_RKMULTIRATE(ts,tab->order,ts->vec_sol,NULL);CHKERRQ(ierr); 788 for (k=0; k<rk->dtratio; k++){ 789 for (i=0; i<s; i++) { 790 rk->stage_time = t + h/rk->dtratio*c[i]; 791 ierr = TSPreStage(ts,rk->stage_time);CHKERRQ(ierr); 792 ierr = VecCopy(prev_sol,Y[i]);CHKERRQ(ierr); 793 ierr = VecCopy(prev_sol,stage_slow);CHKERRQ(ierr); 794 ierr = VecCopy(prev_sol,stage_fast);CHKERRQ(ierr); 795 /* guarantee the stage RHS for slow components do not change while updating the stage RHS for fast components */ 796 ierr = VecCopy(YdotRHS[i],YdotRHS_prev);CHKERRQ(ierr); 797 /* update the fast components stage value */ 798 for (j=0; j<i; j++) w[j] = h/rk->dtratio*A[i*s+j]; 799 ierr = VecMAXPY(stage_fast,i,w,YdotRHS);CHKERRQ(ierr); 800 /* update the slow components stage value */ 801 ierr = VecCopy(prev_sol,Y[0]);CHKERRQ(ierr); 802 ierr = TSInterpolate_RK(ts,t+k*h/rk->dtratio+h/rk->dtratio*c[i],stage_slow);CHKERRQ(ierr); 803 /* combine the update fast components stage value and slow components stage value to Y[i] */ 804 ierr = VecGetArray(Y[i],&y_ptr);CHKERRQ(ierr); 805 ierr = VecGetArray(stage_slow,&slowstage_ptr);CHKERRQ(ierr); 806 ierr = VecGetArray(stage_fast,&faststage_ptr);CHKERRQ(ierr); 807 for (j=0; j<len_slow; j++) y_ptr[is_slow[j]] = slowstage_ptr[is_slow[j]]; 808 for (j=0; j<len_fast; j++) y_ptr[is_fast[j]] = faststage_ptr[is_fast[j]]; 809 ierr = VecRestoreArray(Y[i],&y_ptr);CHKERRQ(ierr); 810 ierr = VecRestoreArray(stage_slow,&slowstage_ptr);CHKERRQ(ierr); 811 ierr = VecRestoreArray(stage_fast,&faststage_ptr);CHKERRQ(ierr); 812 ierr = TSPostStage(ts,rk->stage_time,i,Y); CHKERRQ(ierr); 813 /* compute the stage RHS */ 814 ierr = TSComputeRHSFunction(ts,t+k*h/rk->dtratio+h/rk->dtratio*c[i],Y[i],YdotRHS[i]);CHKERRQ(ierr); 815 /* re-assign the value of slow components in YdotRHS[i] to the calculated value before working on smaller time step-size */ 816 ierr = VecGetArray(YdotRHS[i],&ydotrhs_ptr);CHKERRQ(ierr); 817 ierr = VecGetArray(YdotRHS_prev,&ydotrhsprev_ptr);CHKERRQ(ierr); 818 for (j=0; j<len_slow; j++) ydotrhs_ptr[is_slow[j]] = ydotrhsprev_ptr[is_slow[j]]; 819 ierr = VecRestoreArray(YdotRHS[i],&ydotrhs_ptr);CHKERRQ(ierr); 820 ierr = VecRestoreArray(YdotRHS_prev,&ydotrhsprev_ptr);CHKERRQ(ierr); 821 } 822 rk->slow = 1; 823 ierr = TSEvaluateStep_RKMULTIRATE(ts,tab->order,ts->vec_sol,NULL);CHKERRQ(ierr); 824 /* update the intial value for fast components */ 825 ierr = VecGetSubVector(prev_sol,ts->isf,&presol_fast);CHKERRQ(ierr); 826 ierr = VecGetSubVector(ts->vec_sol,ts->isf,&newslo_fast);CHKERRQ(ierr); 827 ierr = VecCopy(newslo_fast,presol_fast);CHKERRQ(ierr); 828 ierr = VecRestoreSubVector(prev_sol,ts->isf,&presol_fast);CHKERRQ(ierr); 829 ierr = VecRestoreSubVector(ts->vec_sol,ts->isf,&newslo_fast);CHKERRQ(ierr); 830 } 831 ts->ptime += ts->time_step; 832 ts->time_step = next_time_step; 833 rk->slow = 0; 834 rk->status = TS_STEP_COMPLETE; 835 /* free memory of work vectors */ 836 ierr = VecDestroy(&stage_fast);CHKERRQ(ierr); 837 ierr = VecDestroy(&stage_slow);CHKERRQ(ierr); 838 ierr = VecDestroy(&YdotRHS_prev);CHKERRQ(ierr); 839 ierr = VecDestroy(&prev_sol);CHKERRQ(ierr); 840 PetscFunctionReturn(0); 841 } 842 843 /* 844 This is for partitioned RHS multirate RK method 845 The step completion formula is 846 847 x1 = x0 + h b^T YdotRHS 848 849 */ 850 static PetscErrorCode TSEvaluateStep_RKPMULTIRATE(TS ts,PetscInt order,Vec X,PetscBool *done) 851 { 852 TS_RK *rk = (TS_RK*)ts->data; 853 RKTableau tab = rk->tableau; 854 Vec Xslow,Xfast; /* subvectors of X which store slow components and fast components respectively */ 855 PetscScalar *w = rk->work; 856 PetscReal h = ts->time_step; 857 PetscInt s = tab->s,j; 858 PetscErrorCode ierr; 859 860 PetscFunctionBegin; 861 ierr = VecCopy(ts->vec_sol,X);CHKERRQ(ierr); 862 if (!rk->slow){ 863 for (j=0; j<s; j++) w[j] = h*tab->b[j]; 864 ierr = VecGetSubVector(ts->vec_sol,ts->iss,&Xslow);CHKERRQ(ierr); 865 ierr = VecMAXPY(Xslow,s,w,rk->YdotRHS_slow);CHKERRQ(ierr); 866 ierr = VecRestoreSubVector(ts->vec_sol,ts->iss,&Xslow);CHKERRQ(ierr);; 867 } else {; 868 for (j=0; j<s; j++) w[j] = h/rk->dtratio*tab->b[j]; 869 ierr = VecGetSubVector(X,ts->isf,&Xfast);CHKERRQ(ierr); 870 ierr = VecMAXPY(Xfast,s,w,rk->YdotRHS_fast);CHKERRQ(ierr); 871 ierr = VecRestoreSubVector(X,ts->isf,&Xfast);CHKERRQ(ierr); 872 } 873 PetscFunctionReturn(0); 874 } 875 876 static PetscErrorCode TSStep_RKPMULTIRATE(TS ts) 877 { 878 TS_RK *rk = (TS_RK*)ts->data; 879 RKTableau tab = rk->tableau; 880 Vec *Y = rk->Y; 881 Vec *YdotRHS_fast = rk->YdotRHS_fast,*YdotRHS_slow = rk->YdotRHS_slow; 882 Vec Yslow,Yfast; /* subvectors store the stges of slow components and fast components respectively */ 883 Vec Yfast_prev,Yfast_new,prev_sol; /* subvectors store the previous and new solution of fast components and vector store the previous solution */ 884 const PetscInt s = tab->s; 885 const PetscReal *A = tab->A,*c = tab->c; 886 PetscScalar *w = rk->work; 887 PetscInt i,j,k; 888 PetscReal next_time_step = ts->time_step,t = ts->ptime,h = ts->time_step; 889 PetscErrorCode ierr; 890 891 PetscFunctionBegin; 892 rk->status = TS_STEP_INCOMPLETE; 893 for (i=0; i<s; i++) { 894 rk->stage_time = t + h*c[i]; 895 ierr = TSPreStage(ts,rk->stage_time);CHKERRQ(ierr); 896 ierr = VecCopy(ts->vec_sol,Y[i]);CHKERRQ(ierr); 897 /*ierr = VecCopy(ts->vec_sol,Yc[i]);CHKERRQ(ierr);*/ 898 ierr = VecGetSubVector(Y[i],ts->iss,&Yslow);CHKERRQ(ierr); 899 ierr = VecGetSubVector(Y[i],ts->isf,&Yfast);CHKERRQ(ierr); 900 for (j=0; j<i; j++) w[j] = h*A[i*s+j]; 901 ierr = VecMAXPY(Yslow,i,w,YdotRHS_slow);CHKERRQ(ierr); 902 ierr = VecMAXPY(Yfast,i,w,YdotRHS_fast);CHKERRQ(ierr); 903 ierr = VecRestoreSubVector(Y[i],ts->iss,&Yslow);CHKERRQ(ierr); 904 ierr = VecRestoreSubVector(Y[i],ts->isf,&Yfast);CHKERRQ(ierr); 905 ierr = TSPostStage(ts,rk->stage_time,i,Y); CHKERRQ(ierr); 906 /* compute the stage RHS for both slow and fast components */ 907 ierr = TSComputeRHSFunctionslow(ts,t+h*c[i],Y[i],YdotRHS_slow[i]);CHKERRQ(ierr); 908 ierr = TSComputeRHSFunctionfast(ts,t+h*c[i],Y[i],YdotRHS_fast[i]);CHKERRQ(ierr); 909 } 910 /* store the value of slow components at previous time */ 911 ierr = VecDuplicate(ts->vec_sol,&prev_sol);CHKERRQ(ierr); 912 ierr = VecCopy(ts->vec_sol,prev_sol);CHKERRQ(ierr); 913 rk->slow = 0; 914 ierr = TSEvaluateStep_RKPMULTIRATE(ts,tab->order,ts->vec_sol,NULL);CHKERRQ(ierr); 915 for (k=0; k<rk->dtratio; k++){ 916 for (i=0; i<s; i++) { 917 rk->stage_time = t + h/rk->dtratio*c[i]; 918 ierr = TSPreStage(ts,rk->stage_time);CHKERRQ(ierr); 919 ierr = VecCopy(prev_sol,Y[i]);CHKERRQ(ierr); 920 /* stage value for fast components */ 921 for (j=0; j<i; j++) w[j] = h/rk->dtratio*A[i*s+j]; 922 ierr = VecGetSubVector(Y[i],ts->isf,&Yfast);CHKERRQ(ierr); 923 ierr = VecMAXPY(Yfast,i,w,YdotRHS_fast);CHKERRQ(ierr); 924 ierr = VecRestoreSubVector(Y[i],ts->isf,&Yfast);CHKERRQ(ierr); 925 /* stage value for slow components */ 926 ierr = VecCopy(prev_sol,Y[0]);CHKERRQ(ierr); 927 ierr = VecGetSubVector(Y[i],ts->iss,&Yslow);CHKERRQ(ierr); 928 ierr = TSInterpolate_PARTITIONEDMRK(ts,t+k*h/rk->dtratio+h/rk->dtratio*c[i],Yslow);CHKERRQ(ierr); 929 ierr = VecRestoreSubVector(Y[i],ts->iss,&Yslow);CHKERRQ(ierr); 930 ierr = TSPostStage(ts,rk->stage_time,i,Y); CHKERRQ(ierr); 931 /* compute the stage RHS for fast components */ 932 ierr = TSComputeRHSFunctionfast(ts,t+k*h/rk->dtratio+h/rk->dtratio*c[i],Y[i],YdotRHS_fast[i]);CHKERRQ(ierr); 933 } 934 /* update the value of fast components whenusing fast time step */ 935 rk->slow = 1; 936 ierr = TSEvaluateStep_RKPMULTIRATE(ts,tab->order,ts->vec_sol,NULL);CHKERRQ(ierr); 937 ierr = VecGetSubVector(prev_sol,ts->isf,&Yfast_prev);CHKERRQ(ierr); 938 ierr = VecGetSubVector(ts->vec_sol,ts->isf,&Yfast_new);CHKERRQ(ierr); 939 ierr = VecCopy(Yfast_new,Yfast_prev);CHKERRQ(ierr); 940 ierr = VecRestoreSubVector(prev_sol,ts->isf,&Yfast_prev);CHKERRQ(ierr); 941 ierr = VecRestoreSubVector(ts->vec_sol,ts->isf,&Yfast_new);CHKERRQ(ierr); 942 } 943 ts->ptime += ts->time_step; 944 ts->time_step = next_time_step; 945 rk->slow = 0; 946 rk->status = TS_STEP_COMPLETE; 947 ierr = VecDestroy(&prev_sol);CHKERRQ(ierr); 948 PetscFunctionReturn(0); 949 } 950 951 static PetscErrorCode TSAdjointSetUp_RK(TS ts) 952 { 953 TS_RK *rk = (TS_RK*)ts->data; 954 RKTableau tab = rk->tableau; 955 PetscInt s = tab->s; 956 PetscErrorCode ierr; 957 958 PetscFunctionBegin; 959 if (ts->adjointsetupcalled++) PetscFunctionReturn(0); 960 ierr = VecDuplicateVecs(ts->vecs_sensi[0],s*ts->numcost,&rk->VecDeltaLam);CHKERRQ(ierr); 961 if(ts->vecs_sensip) { 962 ierr = VecDuplicateVecs(ts->vecs_sensip[0],s*ts->numcost,&rk->VecDeltaMu);CHKERRQ(ierr); 963 } 964 PetscFunctionReturn(0); 965 } 966 967 static PetscErrorCode TSAdjointStep_RK(TS ts) 968 { 969 TS_RK *rk = (TS_RK*)ts->data; 970 RKTableau tab = rk->tableau; 971 const PetscInt s = tab->s; 972 const PetscReal *A = tab->A,*b = tab->b,*c = tab->c; 973 PetscScalar *w = rk->work; 974 Vec *Y = rk->Y,*VecDeltaLam = rk->VecDeltaLam,*VecDeltaMu = rk->VecDeltaMu; 975 PetscInt i,j,nadj; 976 PetscReal t = ts->ptime; 977 PetscErrorCode ierr; 978 PetscReal h = ts->time_step; 979 980 PetscFunctionBegin; 981 rk->status = TS_STEP_INCOMPLETE; 982 for (i=s-1; i>=0; i--) { 983 Mat J; 984 PetscReal stage_time = t + h*(1.0-c[i]); 985 PetscBool zero = PETSC_FALSE; 986 987 ierr = TSGetRHSJacobian(ts,&J,NULL,NULL,NULL);CHKERRQ(ierr); 988 ierr = TSComputeRHSJacobian(ts,stage_time,Y[i],J,J);CHKERRQ(ierr); 989 if (ts->vec_costintegral) { 990 ierr = TSAdjointComputeDRDYFunction(ts,stage_time,Y[i],ts->vecs_drdy);CHKERRQ(ierr); 991 } 992 /* Stage values of mu */ 993 if (ts->vecs_sensip) { 994 ierr = TSAdjointComputeRHSJacobian(ts,stage_time,Y[i],ts->Jacp);CHKERRQ(ierr); 995 if (ts->vec_costintegral) { 996 ierr = TSAdjointComputeDRDPFunction(ts,stage_time,Y[i],ts->vecs_drdp);CHKERRQ(ierr); 997 } 998 } 999 1000 if (b[i] == 0 && i == s-1) zero = PETSC_TRUE; 1001 for (nadj=0; nadj<ts->numcost; nadj++) { 1002 DM dm; 1003 Vec VecSensiTemp; 1004 1005 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 1006 ierr = DMGetGlobalVector(dm,&VecSensiTemp);CHKERRQ(ierr); 1007 /* Stage values of lambda */ 1008 if (!zero) { 1009 ierr = VecCopy(ts->vecs_sensi[nadj],VecSensiTemp);CHKERRQ(ierr); 1010 ierr = VecScale(VecSensiTemp,-h*b[i]);CHKERRQ(ierr); 1011 for (j=i+1; j<s; j++) w[j-i-1] = -h*A[j*s+i]; 1012 ierr = VecMAXPY(VecSensiTemp,s-i-1,w,&VecDeltaLam[nadj*s+i+1]);CHKERRQ(ierr); 1013 ierr = MatMultTranspose(J,VecSensiTemp,VecDeltaLam[nadj*s+i]);CHKERRQ(ierr); 1014 } else { 1015 ierr = VecSet(VecDeltaLam[nadj*s+i],0);CHKERRQ(ierr); 1016 } 1017 if (ts->vec_costintegral) { 1018 ierr = VecAXPY(VecDeltaLam[nadj*s+i],-h*b[i],ts->vecs_drdy[nadj]);CHKERRQ(ierr); 1019 } 1020 1021 /* Stage values of mu */ 1022 if (ts->vecs_sensip) { 1023 if (!zero) { 1024 ierr = MatMultTranspose(ts->Jacp,VecSensiTemp,VecDeltaMu[nadj*s+i]);CHKERRQ(ierr); 1025 } else { 1026 ierr = VecSet(VecDeltaMu[nadj*s+i],0);CHKERRQ(ierr); 1027 } 1028 if (ts->vec_costintegral) { 1029 ierr = VecAXPY(VecDeltaMu[nadj*s+i],-h*b[i],ts->vecs_drdp[nadj]);CHKERRQ(ierr); 1030 } 1031 } 1032 ierr = DMRestoreGlobalVector(dm,&VecSensiTemp);CHKERRQ(ierr); 1033 } 1034 } 1035 1036 for (j=0; j<s; j++) w[j] = 1.0; 1037 for (nadj=0; nadj<ts->numcost; nadj++) { 1038 ierr = VecMAXPY(ts->vecs_sensi[nadj],s,w,&VecDeltaLam[nadj*s]);CHKERRQ(ierr); 1039 if (ts->vecs_sensip) { 1040 ierr = VecMAXPY(ts->vecs_sensip[nadj],s,w,&VecDeltaMu[nadj*s]);CHKERRQ(ierr); 1041 } 1042 } 1043 rk->status = TS_STEP_COMPLETE; 1044 PetscFunctionReturn(0); 1045 } 1046 1047 static PetscErrorCode TSRKTableauReset(TS ts) 1048 { 1049 TS_RK *rk = (TS_RK*)ts->data; 1050 RKTableau tab = rk->tableau; 1051 PetscErrorCode ierr; 1052 1053 PetscFunctionBegin; 1054 if (!tab) PetscFunctionReturn(0); 1055 ierr = PetscFree(rk->work);CHKERRQ(ierr); 1056 ierr = VecDestroyVecs(tab->s,&rk->Y);CHKERRQ(ierr); 1057 ierr = VecDestroyVecs(tab->s,&rk->YdotRHS);CHKERRQ(ierr); 1058 ierr = VecDestroyVecs(tab->s,&rk->YdotRHS_fast);CHKERRQ(ierr); 1059 ierr = VecDestroyVecs(tab->s,&rk->YdotRHS_slow);CHKERRQ(ierr); 1060 ierr = VecDestroyVecs(tab->s*ts->numcost,&rk->VecDeltaLam);CHKERRQ(ierr); 1061 ierr = VecDestroyVecs(tab->s*ts->numcost,&rk->VecDeltaMu);CHKERRQ(ierr); 1062 PetscFunctionReturn(0); 1063 } 1064 1065 static PetscErrorCode TSReset_RK(TS ts) 1066 { 1067 TS_RK *rk = (TS_RK*)ts->data; 1068 PetscErrorCode ierr; 1069 1070 PetscFunctionBegin; 1071 ierr = TSRKTableauReset(ts);CHKERRQ(ierr); 1072 ierr = VecDestroy(&rk->VecCostIntegral0);CHKERRQ(ierr); 1073 PetscFunctionReturn(0); 1074 } 1075 1076 static PetscErrorCode DMCoarsenHook_TSRK(DM fine,DM coarse,void *ctx) 1077 { 1078 PetscFunctionBegin; 1079 PetscFunctionReturn(0); 1080 } 1081 1082 static PetscErrorCode DMRestrictHook_TSRK(DM fine,Mat restrct,Vec rscale,Mat inject,DM coarse,void *ctx) 1083 { 1084 PetscFunctionBegin; 1085 PetscFunctionReturn(0); 1086 } 1087 1088 1089 static PetscErrorCode DMSubDomainHook_TSRK(DM dm,DM subdm,void *ctx) 1090 { 1091 PetscFunctionBegin; 1092 PetscFunctionReturn(0); 1093 } 1094 1095 static PetscErrorCode DMSubDomainRestrictHook_TSRK(DM dm,VecScatter gscat,VecScatter lscat,DM subdm,void *ctx) 1096 { 1097 1098 PetscFunctionBegin; 1099 PetscFunctionReturn(0); 1100 } 1101 /* 1102 static PetscErrorCode RKSetAdjCoe(RKTableau tab) 1103 { 1104 PetscReal *A,*b; 1105 PetscInt s,i,j; 1106 PetscErrorCode ierr; 1107 1108 PetscFunctionBegin; 1109 s = tab->s; 1110 ierr = PetscMalloc2(s*s,&A,s,&b);CHKERRQ(ierr); 1111 1112 for (i=0; i<s; i++) 1113 for (j=0; j<s; j++) { 1114 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]; 1115 ierr = PetscPrintf(PETSC_COMM_WORLD,"Coefficients: A[%D][%D]=%.6f\n",i,j,A[i*s+j]);CHKERRQ(ierr); 1116 } 1117 1118 for (i=0; i<s; i++) b[i] = (tab->b[s-1-i]==0)? 0: -tab->b[s-1-i]; 1119 1120 ierr = PetscMemcpy(tab->A,A,s*s*sizeof(A[0]));CHKERRQ(ierr); 1121 ierr = PetscMemcpy(tab->b,b,s*sizeof(b[0]));CHKERRQ(ierr); 1122 ierr = PetscFree2(A,b);CHKERRQ(ierr); 1123 PetscFunctionReturn(0); 1124 } 1125 */ 1126 1127 static PetscErrorCode TSRKTableauSetUp(TS ts) 1128 { 1129 TS_RK *rk = (TS_RK*)ts->data; 1130 RKTableau tab = rk->tableau; 1131 Vec YdotRHS_fast,YdotRHS_slow; 1132 PetscErrorCode ierr; 1133 1134 PetscFunctionBegin; 1135 ierr = PetscMalloc1(tab->s,&rk->work);CHKERRQ(ierr); 1136 ierr = VecDuplicateVecs(ts->vec_sol,tab->s,&rk->Y);CHKERRQ(ierr); 1137 ierr = VecDuplicateVecs(ts->vec_sol,tab->s,&rk->YdotRHS);CHKERRQ(ierr); 1138 ierr = VecGetSubVector(ts->vec_sol,ts->isf,&YdotRHS_fast);CHKERRQ(ierr); 1139 ierr = VecDuplicateVecs(YdotRHS_fast,tab->s,&rk->YdotRHS_fast);CHKERRQ(ierr); 1140 ierr = VecRestoreSubVector(ts->vec_sol,ts->isf,&YdotRHS_fast);CHKERRQ(ierr); 1141 ierr = VecGetSubVector(ts->vec_sol,ts->iss,&YdotRHS_slow);CHKERRQ(ierr); 1142 ierr = VecDuplicateVecs(YdotRHS_slow,tab->s,&rk->YdotRHS_slow);CHKERRQ(ierr); 1143 ierr = VecRestoreSubVector(ts->vec_sol,ts->iss,&YdotRHS_slow);CHKERRQ(ierr); 1144 1145 PetscFunctionReturn(0); 1146 } 1147 1148 1149 static PetscErrorCode TSSetUp_RK(TS ts) 1150 { 1151 TS_RK *rk = (TS_RK*)ts->data; 1152 PetscErrorCode ierr; 1153 DM dm; 1154 1155 PetscFunctionBegin; 1156 ierr = TSCheckImplicitTerm(ts);CHKERRQ(ierr); 1157 ierr = TSRKTableauSetUp(ts);CHKERRQ(ierr); 1158 if (!rk->VecCostIntegral0 && ts->vec_costintegral && ts->costintegralfwd) { /* back up cost integral */ 1159 ierr = VecDuplicate(ts->vec_costintegral,&rk->VecCostIntegral0);CHKERRQ(ierr); 1160 } 1161 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 1162 ierr = DMCoarsenHookAdd(dm,DMCoarsenHook_TSRK,DMRestrictHook_TSRK,ts);CHKERRQ(ierr); 1163 ierr = DMSubDomainHookAdd(dm,DMSubDomainHook_TSRK,DMSubDomainRestrictHook_TSRK,ts);CHKERRQ(ierr); 1164 PetscFunctionReturn(0); 1165 } 1166 1167 /* 1168 construct a database to choose single-step rk method, or combined rhs mutirate rk method or partitioned rhs rk method 1169 */ 1170 const char *const TSRKMultirateTypes[] = {"NONE","COMBINED","PARTITIONED","TSRKMultirateType","RKM_",0}; 1171 1172 static PetscErrorCode TSSetFromOptions_RK(PetscOptionItems *PetscOptionsObject,TS ts) 1173 { 1174 TS_RK *rk = (TS_RK*)ts->data; 1175 PetscErrorCode ierr; 1176 1177 PetscFunctionBegin; 1178 ierr = PetscOptionsHead(PetscOptionsObject,"RK ODE solver options");CHKERRQ(ierr); 1179 { 1180 RKTableauLink link; 1181 PetscInt count,choice; 1182 PetscBool flg; 1183 const char **namelist; 1184 PetscInt rkmtype = 0; 1185 1186 for (link=RKTableauList,count=0; link; link=link->next,count++) ; 1187 ierr = PetscMalloc1(count,(char***)&namelist);CHKERRQ(ierr); 1188 for (link=RKTableauList,count=0; link; link=link->next,count++) namelist[count] = link->tab.name; 1189 ierr = PetscOptionsEList("-ts_rk_multirate_type","Use Multirate or partioned RHS Multirate or single rate RK method","TSRKSetMultirateType",TSRKMultirateTypes,3,TSRKMultirateTypes[0],&rkmtype,&flg);CHKERRQ(ierr); 1190 if (flg) {ierr = TSRKSetMultirateType(ts,rkmtype);CHKERRQ(ierr);} 1191 ierr = PetscOptionsEList("-ts_rk_type","Family of RK method","TSRKSetType",(const char*const*)namelist,count,rk->tableau->name,&choice,&flg);CHKERRQ(ierr); 1192 if (flg) {ierr = TSRKSetType(ts,namelist[choice]);CHKERRQ(ierr);} 1193 ierr = PetscFree(namelist);CHKERRQ(ierr); 1194 } 1195 ierr = PetscOptionsTail();CHKERRQ(ierr); 1196 ierr = PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"Multirate methods options","");CHKERRQ(ierr); 1197 ierr = PetscOptionsInt("-ts_rk_dtratio","time step ratio between slow and fast","",rk->dtratio,&rk->dtratio,NULL);CHKERRQ(ierr); 1198 ierr = PetscOptionsEnd();CHKERRQ(ierr); 1199 PetscFunctionReturn(0); 1200 } 1201 1202 static PetscErrorCode TSView_RK(TS ts,PetscViewer viewer) 1203 { 1204 TS_RK *rk = (TS_RK*)ts->data; 1205 PetscBool iascii; 1206 PetscErrorCode ierr; 1207 1208 PetscFunctionBegin; 1209 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 1210 if (iascii) { 1211 RKTableau tab = rk->tableau; 1212 TSRKType rktype; 1213 char buf[512]; 1214 ierr = TSRKGetType(ts,&rktype);CHKERRQ(ierr); 1215 ierr = PetscViewerASCIIPrintf(viewer," RK type %s\n",rktype);CHKERRQ(ierr); 1216 ierr = PetscViewerASCIIPrintf(viewer," Order: %D\n",tab->order);CHKERRQ(ierr); 1217 ierr = PetscViewerASCIIPrintf(viewer," FSAL property: %s\n",tab->FSAL ? "yes" : "no");CHKERRQ(ierr); 1218 ierr = PetscFormatRealArray(buf,sizeof(buf),"% 8.6f",tab->s,tab->c);CHKERRQ(ierr); 1219 ierr = PetscViewerASCIIPrintf(viewer," Abscissa c = %s\n",buf);CHKERRQ(ierr); 1220 } 1221 PetscFunctionReturn(0); 1222 } 1223 1224 static PetscErrorCode TSLoad_RK(TS ts,PetscViewer viewer) 1225 { 1226 PetscErrorCode ierr; 1227 TSAdapt adapt; 1228 1229 PetscFunctionBegin; 1230 ierr = TSGetAdapt(ts,&adapt);CHKERRQ(ierr); 1231 ierr = TSAdaptLoad(adapt,viewer);CHKERRQ(ierr); 1232 PetscFunctionReturn(0); 1233 } 1234 1235 /*@C 1236 TSRKSetType - Set the type of RK scheme 1237 1238 Logically collective 1239 1240 Input Parameter: 1241 + ts - timestepping context 1242 - rktype - type of RK-scheme 1243 1244 Options Database: 1245 . -ts_rk_type - <1fe,2a,3,3bs,4,5f,5dp,5bs> 1246 1247 Level: intermediate 1248 1249 .seealso: TSRKGetType(), TSRK, TSRKType, TSRK1FE, TSRK2A, TSRK3, TSRK3BS, TSRK4, TSRK5F, TSRK5DP, TSRK5BS 1250 @*/ 1251 PetscErrorCode TSRKSetType(TS ts,TSRKType rktype) 1252 { 1253 PetscErrorCode ierr; 1254 1255 PetscFunctionBegin; 1256 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1257 PetscValidCharPointer(rktype,2); 1258 ierr = PetscTryMethod(ts,"TSRKSetType_C",(TS,TSRKType),(ts,rktype));CHKERRQ(ierr); 1259 PetscFunctionReturn(0); 1260 } 1261 1262 /*@C 1263 TSRKSetMultirateType - Set the type of RK Multirate scheme 1264 1265 Logically collective 1266 1267 Input Parameter: 1268 + ts - timestepping context 1269 - rkmtype - type of RKM-scheme 1270 1271 Options Database: 1272 . -ts_rk_multiarte_type - <none,combined,partitioned> 1273 1274 Level: intermediate 1275 @*/ 1276 PetscErrorCode TSRKSetMultirateType(TS ts, TSRKMultirateType rkmtype) 1277 { 1278 TS_RK *rk = (TS_RK*)ts->data; 1279 PetscErrorCode ierr; 1280 1281 PetscFunctionBegin; 1282 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1283 switch(rkmtype){ 1284 case RKM_NONE: 1285 break; 1286 case RKM_COMBINED: 1287 ts->ops->step = TSStep_RKMULTIRATE; 1288 ts->ops->evaluatestep = TSEvaluateStep_RKMULTIRATE; 1289 rk->dtratio = 2; 1290 ierr = TSRKSetType(ts,TSMRKDefault);CHKERRQ(ierr); 1291 break; 1292 case RKM_PARTITIONED: 1293 ts->ops->step = TSStep_RKPMULTIRATE; 1294 ts->ops->evaluatestep = TSEvaluateStep_RKPMULTIRATE; 1295 ts->ops->interpolate = TSInterpolate_PARTITIONEDMRK; 1296 rk->dtratio = 2; 1297 ierr = TSRKSetType(ts,TSMRKDefault);CHKERRQ(ierr); 1298 break; 1299 default : 1300 SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown type '%s'",rkmtype); 1301 } 1302 PetscFunctionReturn(0); 1303 } 1304 1305 /*@C 1306 TSRKGetType - Get the type of RK scheme 1307 1308 Logically collective 1309 1310 Input Parameter: 1311 . ts - timestepping context 1312 1313 Output Parameter: 1314 . rktype - type of RK-scheme 1315 1316 Level: intermediate 1317 1318 .seealso: TSRKGetType() 1319 @*/ 1320 PetscErrorCode TSRKGetType(TS ts,TSRKType *rktype) 1321 { 1322 PetscErrorCode ierr; 1323 1324 PetscFunctionBegin; 1325 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1326 ierr = PetscUseMethod(ts,"TSRKGetType_C",(TS,TSRKType*),(ts,rktype));CHKERRQ(ierr); 1327 PetscFunctionReturn(0); 1328 } 1329 1330 static PetscErrorCode TSRKGetType_RK(TS ts,TSRKType *rktype) 1331 { 1332 TS_RK *rk = (TS_RK*)ts->data; 1333 1334 PetscFunctionBegin; 1335 *rktype = rk->tableau->name; 1336 PetscFunctionReturn(0); 1337 } 1338 1339 static PetscErrorCode TSRKSetType_RK(TS ts,TSRKType rktype) 1340 { 1341 TS_RK *rk = (TS_RK*)ts->data; 1342 PetscErrorCode ierr; 1343 PetscBool match; 1344 RKTableauLink link; 1345 1346 PetscFunctionBegin; 1347 if (rk->tableau) { 1348 ierr = PetscStrcmp(rk->tableau->name,rktype,&match);CHKERRQ(ierr); 1349 if (match) PetscFunctionReturn(0); 1350 } 1351 for (link = RKTableauList; link; link=link->next) { 1352 ierr = PetscStrcmp(link->tab.name,rktype,&match);CHKERRQ(ierr); 1353 if (match) { 1354 if (ts->setupcalled) {ierr = TSRKTableauReset(ts);CHKERRQ(ierr);} 1355 rk->tableau = &link->tab; 1356 if (ts->setupcalled) {ierr = TSRKTableauSetUp(ts);CHKERRQ(ierr);} 1357 ts->default_adapt_type = rk->tableau->bembed ? TSADAPTBASIC : TSADAPTNONE; 1358 PetscFunctionReturn(0); 1359 } 1360 } 1361 SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_UNKNOWN_TYPE,"Could not find '%s'",rktype); 1362 PetscFunctionReturn(0); 1363 } 1364 1365 static PetscErrorCode TSGetStages_RK(TS ts,PetscInt *ns,Vec **Y) 1366 { 1367 TS_RK *rk = (TS_RK*)ts->data; 1368 1369 PetscFunctionBegin; 1370 if (ns) *ns = rk->tableau->s; 1371 if (Y) *Y = rk->Y; 1372 PetscFunctionReturn(0); 1373 } 1374 1375 static PetscErrorCode TSDestroy_RK(TS ts) 1376 { 1377 PetscErrorCode ierr; 1378 1379 PetscFunctionBegin; 1380 ierr = TSReset_RK(ts);CHKERRQ(ierr); 1381 if (ts->dm) { 1382 ierr = DMCoarsenHookRemove(ts->dm,DMCoarsenHook_TSRK,DMRestrictHook_TSRK,ts);CHKERRQ(ierr); 1383 ierr = DMSubDomainHookRemove(ts->dm,DMSubDomainHook_TSRK,DMSubDomainRestrictHook_TSRK,ts);CHKERRQ(ierr); 1384 } 1385 ierr = PetscFree(ts->data);CHKERRQ(ierr); 1386 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSRKGetType_C",NULL);CHKERRQ(ierr); 1387 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSRKSetType_C",NULL);CHKERRQ(ierr); 1388 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSRKSetMultirateType_C",NULL);CHKERRQ(ierr); 1389 PetscFunctionReturn(0); 1390 } 1391 1392 /*MC 1393 TSRK - ODE and DAE solver using Runge-Kutta schemes 1394 1395 The user should provide the right hand side of the equation 1396 using TSSetRHSFunction(). 1397 1398 Notes: 1399 The default is TSRK3BS, it can be changed with TSRKSetType() or -ts_rk_type 1400 1401 Level: beginner 1402 1403 .seealso: TSCreate(), TS, TSSetType(), TSRKSetType(), TSRKGetType(), TSRKSetFullyImplicit(), TSRK2D, TTSRK2E, TSRK3, 1404 TSRK4, TSRK5, TSRKPRSSP2, TSRKBPR3, TSRKType, TSRKRegister(), TSRKSetMultirateType() 1405 1406 M*/ 1407 PETSC_EXTERN PetscErrorCode TSCreate_RK(TS ts) 1408 { 1409 TS_RK *rk; 1410 PetscErrorCode ierr; 1411 1412 PetscFunctionBegin; 1413 ierr = TSRKInitializePackage();CHKERRQ(ierr); 1414 1415 ts->ops->reset = TSReset_RK; 1416 ts->ops->destroy = TSDestroy_RK; 1417 ts->ops->view = TSView_RK; 1418 ts->ops->load = TSLoad_RK; 1419 ts->ops->setup = TSSetUp_RK; 1420 ts->ops->adjointsetup = TSAdjointSetUp_RK; 1421 ts->ops->interpolate = TSInterpolate_RK; 1422 ts->ops->step = TSStep_RK; 1423 ts->ops->evaluatestep = TSEvaluateStep_RK; 1424 ts->ops->rollback = TSRollBack_RK; 1425 ts->ops->setfromoptions = TSSetFromOptions_RK; 1426 ts->ops->getstages = TSGetStages_RK; 1427 ts->ops->adjointstep = TSAdjointStep_RK; 1428 1429 ts->ops->adjointintegral = TSAdjointCostIntegral_RK; 1430 ts->ops->forwardintegral = TSForwardCostIntegral_RK; 1431 1432 ierr = PetscNewLog(ts,&rk);CHKERRQ(ierr); 1433 ts->data = (void*)rk; 1434 1435 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSRKGetType_C",TSRKGetType_RK);CHKERRQ(ierr); 1436 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSRKSetType_C",TSRKSetType_RK);CHKERRQ(ierr); 1437 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSRKSetMultirateType_C",TSRKSetMultirateType);CHKERRQ(ierr); 1438 1439 ierr = TSRKSetType(ts,TSRKDefault);CHKERRQ(ierr); 1440 PetscFunctionReturn(0); 1441 } 1442