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