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 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)}, 250 {0,0,0,0,0}, 251 {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)}, 252 {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)}, 253 {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)}, 254 {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)}, 255 {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)}}; 256 257 ierr = TSRKRegister(TSRK5DP,5,7,&A[0][0],b,NULL,bembed,5,binterp[0]);CHKERRQ(ierr); 258 } 259 { 260 const PetscReal 261 A[8][8] = {{0,0,0,0,0,0,0,0}, 262 {RC(1.0)/RC(6.0),0,0,0,0,0,0,0}, 263 {RC(2.0)/RC(27.0),RC(4.0)/RC(27.0),0,0,0,0,0,0}, 264 {RC(183.0)/RC(1372.0),RC(-162.0)/RC(343.0),RC(1053.0)/RC(1372.0),0,0,0,0,0}, 265 {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}, 266 {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}, 267 {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}, 268 {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}}, 269 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}, 270 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)}; 271 ierr = TSRKRegister(TSRK5BS,5,8,&A[0][0],b,NULL,bembed,0,NULL);CHKERRQ(ierr); 272 } 273 #undef RC 274 PetscFunctionReturn(0); 275 } 276 277 /*@C 278 TSRKRegisterDestroy - Frees the list of schemes that were registered by TSRKRegister(). 279 280 Not Collective 281 282 Level: advanced 283 284 .keywords: TSRK, register, destroy 285 .seealso: TSRKRegister(), TSRKRegisterAll() 286 @*/ 287 PetscErrorCode TSRKRegisterDestroy(void) 288 { 289 PetscErrorCode ierr; 290 RKTableauLink link; 291 292 PetscFunctionBegin; 293 while ((link = RKTableauList)) { 294 RKTableau t = &link->tab; 295 RKTableauList = link->next; 296 ierr = PetscFree3(t->A,t->b,t->c); CHKERRQ(ierr); 297 ierr = PetscFree (t->bembed); CHKERRQ(ierr); 298 ierr = PetscFree (t->binterp); CHKERRQ(ierr); 299 ierr = PetscFree (t->name); CHKERRQ(ierr); 300 ierr = PetscFree (link); CHKERRQ(ierr); 301 } 302 TSRKRegisterAllCalled = PETSC_FALSE; 303 PetscFunctionReturn(0); 304 } 305 306 /*@C 307 TSRKInitializePackage - This function initializes everything in the TSRK package. It is called 308 from TSInitializePackage(). 309 310 Level: developer 311 312 .keywords: TS, TSRK, initialize, package 313 .seealso: PetscInitialize() 314 @*/ 315 PetscErrorCode TSRKInitializePackage(void) 316 { 317 PetscErrorCode ierr; 318 319 PetscFunctionBegin; 320 if (TSRKPackageInitialized) PetscFunctionReturn(0); 321 TSRKPackageInitialized = PETSC_TRUE; 322 ierr = TSRKRegisterAll();CHKERRQ(ierr); 323 ierr = PetscRegisterFinalize(TSRKFinalizePackage);CHKERRQ(ierr); 324 PetscFunctionReturn(0); 325 } 326 327 /*@C 328 TSRKFinalizePackage - This function destroys everything in the TSRK package. It is 329 called from PetscFinalize(). 330 331 Level: developer 332 333 .keywords: Petsc, destroy, package 334 .seealso: PetscFinalize() 335 @*/ 336 PetscErrorCode TSRKFinalizePackage(void) 337 { 338 PetscErrorCode ierr; 339 340 PetscFunctionBegin; 341 TSRKPackageInitialized = PETSC_FALSE; 342 ierr = TSRKRegisterDestroy();CHKERRQ(ierr); 343 PetscFunctionReturn(0); 344 } 345 346 /*@C 347 TSRKRegister - register an RK scheme by providing the entries in the Butcher tableau and optionally embedded approximations and interpolation 348 349 Not Collective, but the same schemes should be registered on all processes on which they will be used 350 351 Input Parameters: 352 + name - identifier for method 353 . order - approximation order of method 354 . s - number of stages, this is the dimension of the matrices below 355 . A - stage coefficients (dimension s*s, row-major) 356 . b - step completion table (dimension s; NULL to use last row of A) 357 . c - abscissa (dimension s; NULL to use row sums of A) 358 . bembed - completion table for embedded method (dimension s; NULL if not available) 359 . p - Order of the interpolation scheme, equal to the number of columns of binterp 360 - binterp - Coefficients of the interpolation formula (dimension s*p; NULL to reuse b with p=1) 361 362 Notes: 363 Several RK methods are provided, this function is only needed to create new methods. 364 365 Level: advanced 366 367 .keywords: TS, register 368 369 .seealso: TSRK 370 @*/ 371 PetscErrorCode TSRKRegister(TSRKType name,PetscInt order,PetscInt s, 372 const PetscReal A[],const PetscReal b[],const PetscReal c[], 373 const PetscReal bembed[],PetscInt p,const PetscReal binterp[]) 374 { 375 PetscErrorCode ierr; 376 RKTableauLink link; 377 RKTableau t; 378 PetscInt i,j; 379 380 PetscFunctionBegin; 381 PetscValidCharPointer(name,1); 382 PetscValidRealPointer(A,4); 383 if (b) PetscValidRealPointer(b,5); 384 if (c) PetscValidRealPointer(c,6); 385 if (bembed) PetscValidRealPointer(bembed,7); 386 if (binterp || p > 1) PetscValidRealPointer(binterp,9); 387 388 ierr = TSRKInitializePackage();CHKERRQ(ierr); 389 ierr = PetscNew(&link);CHKERRQ(ierr); 390 t = &link->tab; 391 392 ierr = PetscStrallocpy(name,&t->name);CHKERRQ(ierr); 393 t->order = order; 394 t->s = s; 395 ierr = PetscMalloc3(s*s,&t->A,s,&t->b,s,&t->c);CHKERRQ(ierr); 396 ierr = PetscMemcpy(t->A,A,s*s*sizeof(A[0]));CHKERRQ(ierr); 397 if (b) { ierr = PetscMemcpy(t->b,b,s*sizeof(b[0]));CHKERRQ(ierr); } 398 else for (i=0; i<s; i++) t->b[i] = A[(s-1)*s+i]; 399 if (c) { ierr = PetscMemcpy(t->c,c,s*sizeof(c[0]));CHKERRQ(ierr); } 400 else for (i=0; i<s; i++) for (j=0,t->c[i]=0; j<s; j++) t->c[i] += A[i*s+j]; 401 t->FSAL = PETSC_TRUE; 402 for (i=0; i<s; i++) if (t->A[(s-1)*s+i] != t->b[i]) t->FSAL = PETSC_FALSE; 403 404 if (bembed) { 405 ierr = PetscMalloc1(s,&t->bembed);CHKERRQ(ierr); 406 ierr = PetscMemcpy(t->bembed,bembed,s*sizeof(bembed[0]));CHKERRQ(ierr); 407 } 408 409 if (!binterp) { p = 1; binterp = t->b; } 410 t->p = p; 411 ierr = PetscMalloc1(s*p,&t->binterp);CHKERRQ(ierr); 412 ierr = PetscMemcpy(t->binterp,binterp,s*p*sizeof(binterp[0]));CHKERRQ(ierr); 413 414 link->next = RKTableauList; 415 RKTableauList = link; 416 PetscFunctionReturn(0); 417 } 418 419 /* 420 This is for single-step RK method 421 The step completion formula is 422 423 x1 = x0 + h b^T YdotRHS 424 425 This function can be called before or after ts->vec_sol has been updated. 426 Suppose we have a completion formula (b) and an embedded formula (be) of different order. 427 We can write 428 429 x1e = x0 + h be^T YdotRHS 430 = x1 - h b^T YdotRHS + h be^T YdotRHS 431 = x1 + h (be - b)^T YdotRHS 432 433 so we can evaluate the method with different order even after the step has been optimistically completed. 434 */ 435 static PetscErrorCode TSEvaluateStep_RK(TS ts,PetscInt order,Vec X,PetscBool *done) 436 { 437 TS_RK *rk = (TS_RK*)ts->data; 438 RKTableau tab = rk->tableau; 439 PetscScalar *w = rk->work; 440 PetscReal h; 441 PetscInt s = tab->s,j; 442 PetscErrorCode ierr; 443 444 PetscFunctionBegin; 445 switch (rk->status) { 446 case TS_STEP_INCOMPLETE: 447 case TS_STEP_PENDING: 448 h = ts->time_step; break; 449 case TS_STEP_COMPLETE: 450 h = ts->ptime - ts->ptime_prev; break; 451 default: SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_PLIB,"Invalid TSStepStatus"); 452 } 453 if (order == tab->order) { 454 if (rk->status == TS_STEP_INCOMPLETE) { 455 ierr = VecCopy(ts->vec_sol,X);CHKERRQ(ierr); 456 for (j=0; j<s; j++) w[j] = h*tab->b[j]; 457 ierr = VecMAXPY(X,s,w,rk->YdotRHS);CHKERRQ(ierr); 458 } else {ierr = VecCopy(ts->vec_sol,X);CHKERRQ(ierr);} 459 PetscFunctionReturn(0); 460 } else if (order == tab->order-1) { 461 if (!tab->bembed) goto unavailable; 462 if (rk->status == TS_STEP_INCOMPLETE) { /*Complete with the embedded method (be)*/ 463 ierr = VecCopy(ts->vec_sol,X);CHKERRQ(ierr); 464 for (j=0; j<s; j++) w[j] = h*tab->bembed[j]; 465 ierr = VecMAXPY(X,s,w,rk->YdotRHS);CHKERRQ(ierr); 466 } else { /*Rollback and re-complete using (be-b) */ 467 ierr = VecCopy(ts->vec_sol,X);CHKERRQ(ierr); 468 for (j=0; j<s; j++) w[j] = h*(tab->bembed[j] - tab->b[j]); 469 ierr = VecMAXPY(X,s,w,rk->YdotRHS);CHKERRQ(ierr); 470 } 471 if (done) *done = PETSC_TRUE; 472 PetscFunctionReturn(0); 473 } 474 unavailable: 475 if (done) *done = PETSC_FALSE; 476 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); 477 PetscFunctionReturn(0); 478 } 479 480 static PetscErrorCode TSForwardCostIntegral_RK(TS ts) 481 { 482 TS_RK *rk = (TS_RK*)ts->data; 483 RKTableau tab = rk->tableau; 484 const PetscInt s = tab->s; 485 const PetscReal *b = tab->b,*c = tab->c; 486 Vec *Y = rk->Y; 487 PetscInt i; 488 PetscErrorCode ierr; 489 490 PetscFunctionBegin; 491 /* backup cost integral */ 492 ierr = VecCopy(ts->vec_costintegral,rk->VecCostIntegral0);CHKERRQ(ierr); 493 for (i=s-1; i>=0; i--) { 494 /* Evolve ts->vec_costintegral to compute integrals */ 495 ierr = TSComputeCostIntegrand(ts,rk->ptime+rk->time_step*c[i],Y[i],ts->vec_costintegrand);CHKERRQ(ierr); 496 ierr = VecAXPY(ts->vec_costintegral,rk->time_step*b[i],ts->vec_costintegrand);CHKERRQ(ierr); 497 } 498 PetscFunctionReturn(0); 499 } 500 501 static PetscErrorCode TSAdjointCostIntegral_RK(TS ts) 502 { 503 TS_RK *rk = (TS_RK*)ts->data; 504 RKTableau tab = rk->tableau; 505 const PetscInt s = tab->s; 506 const PetscReal *b = tab->b,*c = tab->c; 507 Vec *Y = rk->Y; 508 PetscInt i; 509 PetscErrorCode ierr; 510 511 PetscFunctionBegin; 512 for (i=s-1; i>=0; i--) { 513 /* Evolve ts->vec_costintegral to compute integrals */ 514 ierr = TSComputeCostIntegrand(ts,ts->ptime+ts->time_step*(1.0-c[i]),Y[i],ts->vec_costintegrand);CHKERRQ(ierr); 515 ierr = VecAXPY(ts->vec_costintegral,-ts->time_step*b[i],ts->vec_costintegrand);CHKERRQ(ierr); 516 } 517 PetscFunctionReturn(0); 518 } 519 520 static PetscErrorCode TSRollBack_RK(TS ts) 521 { 522 TS_RK *rk = (TS_RK*)ts->data; 523 RKTableau tab = rk->tableau; 524 const PetscInt s = tab->s; 525 const PetscReal *b = tab->b; 526 PetscScalar *w = rk->work; 527 Vec *YdotRHS = rk->YdotRHS; 528 PetscInt j; 529 PetscReal h; 530 PetscErrorCode ierr; 531 532 PetscFunctionBegin; 533 switch (rk->status) { 534 case TS_STEP_INCOMPLETE: 535 case TS_STEP_PENDING: 536 h = ts->time_step; break; 537 case TS_STEP_COMPLETE: 538 h = ts->ptime - ts->ptime_prev; break; 539 default: SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_PLIB,"Invalid TSStepStatus"); 540 } 541 for (j=0; j<s; j++) w[j] = -h*b[j]; 542 ierr = VecMAXPY(ts->vec_sol,s,w,YdotRHS);CHKERRQ(ierr); 543 PetscFunctionReturn(0); 544 } 545 546 static PetscErrorCode TSStep_RK(TS ts) 547 { 548 TS_RK *rk = (TS_RK*)ts->data; 549 RKTableau tab = rk->tableau; 550 const PetscInt s = tab->s; 551 const PetscReal *A = tab->A,*c = tab->c; 552 PetscScalar *w = rk->work; 553 Vec *Y = rk->Y,*YdotRHS = rk->YdotRHS; 554 PetscBool FSAL = tab->FSAL; 555 TSAdapt adapt; 556 PetscInt i,j; 557 PetscInt rejections = 0; 558 PetscBool stageok,accept = PETSC_TRUE; 559 PetscReal next_time_step = ts->time_step; 560 PetscErrorCode ierr; 561 562 PetscFunctionBegin; 563 if (ts->steprollback || ts->steprestart) FSAL = PETSC_FALSE; 564 if (FSAL) { ierr = VecCopy(YdotRHS[s-1],YdotRHS[0]);CHKERRQ(ierr); } 565 566 rk->status = TS_STEP_INCOMPLETE; 567 while (!ts->reason && rk->status != TS_STEP_COMPLETE) { 568 PetscReal t = ts->ptime; 569 PetscReal h = ts->time_step; 570 for (i=0; i<s; i++) { 571 rk->stage_time = t + h*c[i]; 572 ierr = TSPreStage(ts,rk->stage_time); CHKERRQ(ierr); 573 ierr = VecCopy(ts->vec_sol,Y[i]);CHKERRQ(ierr); 574 for (j=0; j<i; j++) w[j] = h*A[i*s+j]; 575 ierr = VecMAXPY(Y[i],i,w,YdotRHS);CHKERRQ(ierr); 576 ierr = TSPostStage(ts,rk->stage_time,i,Y); CHKERRQ(ierr); 577 ierr = TSGetAdapt(ts,&adapt);CHKERRQ(ierr); 578 ierr = TSAdaptCheckStage(adapt,ts,rk->stage_time,Y[i],&stageok);CHKERRQ(ierr); 579 if (!stageok) goto reject_step; 580 if (FSAL && !i) continue; 581 ierr = TSComputeRHSFunction(ts,t+h*c[i],Y[i],YdotRHS[i]);CHKERRQ(ierr); 582 } 583 584 rk->status = TS_STEP_INCOMPLETE; 585 ierr = TSEvaluateStep(ts,tab->order,ts->vec_sol,NULL);CHKERRQ(ierr); 586 rk->status = TS_STEP_PENDING; 587 ierr = TSGetAdapt(ts,&adapt);CHKERRQ(ierr); 588 ierr = TSAdaptCandidatesClear(adapt);CHKERRQ(ierr); 589 ierr = TSAdaptCandidateAdd(adapt,tab->name,tab->order,1,tab->ccfl,(PetscReal)tab->s,PETSC_TRUE);CHKERRQ(ierr); 590 ierr = TSAdaptChoose(adapt,ts,ts->time_step,NULL,&next_time_step,&accept);CHKERRQ(ierr); 591 rk->status = accept ? TS_STEP_COMPLETE : TS_STEP_INCOMPLETE; 592 if (!accept) { /*Roll back the current step */ 593 ierr = TSRollBack_RK(ts);CHKERRQ(ierr); 594 ts->time_step = next_time_step; 595 goto reject_step; 596 } 597 598 if (ts->costintegralfwd) { /*Save the info for the later use in cost integral evaluation*/ 599 rk->ptime = ts->ptime; 600 rk->time_step = ts->time_step; 601 } 602 603 ts->ptime += ts->time_step; 604 ts->time_step = next_time_step; 605 break; 606 607 reject_step: 608 ts->reject++; accept = PETSC_FALSE; 609 if (!ts->reason && ++rejections > ts->max_reject && ts->max_reject >= 0) { 610 ts->reason = TS_DIVERGED_STEP_REJECTED; 611 ierr = PetscInfo2(ts,"Step=%D, step rejections %D greater than current TS allowed, stopping solve\n",ts->steps,rejections);CHKERRQ(ierr); 612 } 613 } 614 PetscFunctionReturn(0); 615 } 616 617 static PetscErrorCode TSInterpolate_RK(TS ts,PetscReal itime,Vec X) 618 { 619 TS_RK *rk = (TS_RK*)ts->data; 620 PetscInt s = rk->tableau->s,p = rk->tableau->p,i,j; 621 PetscReal h; 622 PetscReal tt,t; 623 PetscScalar *b; 624 const PetscReal *B = rk->tableau->binterp; 625 PetscErrorCode ierr; 626 627 PetscFunctionBegin; 628 if (!B) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"TSRK %s does not have an interpolation formula",rk->tableau->name); 629 630 switch (rk->status) { 631 case TS_STEP_INCOMPLETE: 632 case TS_STEP_PENDING: 633 h = ts->time_step; 634 t = (itime - ts->ptime)/h; 635 break; 636 case TS_STEP_COMPLETE: 637 h = ts->ptime - ts->ptime_prev; 638 t = (itime - ts->ptime)/h + 1; /* In the interval [0,1] */ 639 break; 640 default: SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_PLIB,"Invalid TSStepStatus"); 641 } 642 ierr = PetscMalloc1(s,&b);CHKERRQ(ierr); 643 for (i=0; i<s; i++) b[i] = 0; 644 for (j=0,tt=t; j<p; j++,tt*=t) { 645 for (i=0; i<s; i++) { 646 b[i] += h * B[i*p+j] * tt; 647 } 648 } 649 ierr = VecCopy(rk->Y[0],X);CHKERRQ(ierr); 650 ierr = VecMAXPY(X,s,b,rk->YdotRHS);CHKERRQ(ierr); 651 ierr = PetscFree(b);CHKERRQ(ierr); 652 PetscFunctionReturn(0); 653 } 654 655 /* 656 This is interpolate formula for partitioned RHS multirate RK method 657 */ 658 659 static PetscErrorCode TSInterpolate_PARTITIONEDMRK(TS ts,PetscReal itime,Vec X) 660 { 661 TS_RK *rk = (TS_RK*)ts->data; 662 PetscInt s = rk->tableau->s,p = rk->tableau->p,i,j; 663 Vec Yslow; /* vector holds the slow components in Y[0] */ 664 PetscReal h; 665 PetscReal tt,t; 666 PetscScalar *b; 667 const PetscReal *B = rk->tableau->binterp; 668 PetscErrorCode ierr; 669 670 PetscFunctionBegin; 671 if (!B) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"TSRK %s does not have an interpolation formula",rk->tableau->name); 672 673 switch (rk->status) { 674 case TS_STEP_INCOMPLETE: 675 case TS_STEP_PENDING: 676 h = ts->time_step; 677 t = (itime - ts->ptime)/h; 678 break; 679 case TS_STEP_COMPLETE: 680 h = ts->ptime - ts->ptime_prev; 681 t = (itime - ts->ptime)/h + 1; /* In the interval [0,1] */ 682 break; 683 default: SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_PLIB,"Invalid TSStepStatus"); 684 } 685 ierr = PetscMalloc1(s,&b);CHKERRQ(ierr); 686 for (i=0; i<s; i++) b[i] = 0; 687 for (j=0,tt=t; j<p; j++,tt*=t) { 688 for (i=0; i<s; i++) { 689 b[i] += h * B[i*p+j] * tt; 690 } 691 } 692 ierr = VecGetSubVector(rk->Y[0],ts->iss,&Yslow);CHKERRQ(ierr); 693 ierr = VecCopy(Yslow,X);CHKERRQ(ierr); 694 ierr = VecMAXPY(X,s,b,rk->YdotRHS_slow);CHKERRQ(ierr); 695 ierr = VecRestoreSubVector(rk->Y[0],ts->iss,&Yslow);CHKERRQ(ierr); 696 ierr = PetscFree(b);CHKERRQ(ierr); 697 PetscFunctionReturn(0); 698 } 699 700 /* 701 This is for combined RHS multirate RK method 702 The step completion formula is 703 704 x1 = x0 + h b^T YdotRHS 705 706 */ 707 static PetscErrorCode TSEvaluateStep_RKMULTIRATE(TS ts,PetscInt order,Vec X,PetscBool *done) 708 { 709 TS_RK *rk = (TS_RK*)ts->data; 710 RKTableau tab = rk->tableau; 711 Vec Xtemp; /* temporary work vector for X */ 712 PetscScalar *w = rk->work; 713 PetscScalar *x_ptr,*xtemp_ptr; /* location to put the pointer to X and Xtemp respectively */ 714 PetscReal h = ts->time_step; 715 PetscInt s = tab->s,j; 716 PetscInt len_slow,len_fast; /* the number of slow components and fast components respectively */ 717 const PetscInt *is_slow,*is_fast; /* the index for slow components and fast components respectively */ 718 PetscErrorCode ierr; 719 720 PetscFunctionBegin; 721 ierr = VecCopy(ts->vec_sol,X);CHKERRQ(ierr); 722 ierr = VecDuplicate(ts->vec_sol,&Xtemp);CHKERRQ(ierr); 723 ierr = VecCopy(ts->vec_sol,Xtemp);CHKERRQ(ierr); 724 if (!rk->slow){ 725 for (j=0; j<s; j++) w[j] = h*tab->b[j]; 726 /* update the value of slow components, and discard the updating value of fast components */ 727 ierr = VecMAXPY(Xtemp,s,w,rk->YdotRHS);CHKERRQ(ierr); 728 ierr = VecGetArray(X,&x_ptr);CHKERRQ(ierr); 729 ierr = VecGetArray(Xtemp,&xtemp_ptr);CHKERRQ(ierr); 730 ierr = ISGetSize(ts->iss,&len_slow);CHKERRQ(ierr); 731 ierr = ISGetIndices(ts->iss,&is_slow);CHKERRQ(ierr); 732 ierr = ISRestoreIndices(ts->iss,&is_slow);CHKERRQ(ierr); 733 for (j=0; j<len_slow; j++) x_ptr[is_slow[j]] = xtemp_ptr[is_slow[j]]; 734 ierr = VecRestoreArray(X,&x_ptr);CHKERRQ(ierr); 735 ierr = VecRestoreArray(Xtemp,&xtemp_ptr);CHKERRQ(ierr); 736 } else { 737 /* update the value of fast components, and discard the updating value of slow components */ 738 for (j=0; j<s; j++) w[j] = h/rk->dtratio*tab->b[j]; 739 ierr = VecMAXPY(Xtemp,s,w,rk->YdotRHS);CHKERRQ(ierr); 740 ierr = VecGetArray(X,&x_ptr);CHKERRQ(ierr); 741 ierr = VecGetArray(Xtemp,&xtemp_ptr);CHKERRQ(ierr); 742 ierr = ISGetSize(ts->isf,&len_fast);CHKERRQ(ierr); 743 ierr = ISGetIndices(ts->isf,&is_fast);CHKERRQ(ierr); 744 ierr = ISRestoreIndices(ts->isf,&is_fast);CHKERRQ(ierr); 745 for (j=0; j<len_fast; j++) x_ptr[is_fast[j]] = xtemp_ptr[is_fast[j]]; 746 ierr = VecRestoreArray(X,&x_ptr);CHKERRQ(ierr); 747 ierr = VecRestoreArray(Xtemp,&xtemp_ptr);CHKERRQ(ierr); 748 } 749 /* free temporary vector space */ 750 ierr = VecDestroy(&Xtemp);CHKERRQ(ierr); 751 PetscFunctionReturn(0); 752 } 753 754 static PetscErrorCode TSStep_RKMULTIRATE(TS ts) 755 { 756 TS_RK *rk = (TS_RK*)ts->data; 757 RKTableau tab = rk->tableau; 758 Vec *Y = rk->Y,*YdotRHS = rk->YdotRHS; 759 Vec stage_fast,stage_slow; /* vectors store the stage value of fast and slow components respectively */ 760 Vec presol_fast,newslo_fast; /* vectors store the previous and new time solution of fast components respectively */ 761 Vec YdotRHS_prev,prev_sol; /* vectors store the previous YdotRHS and solution respectively */ 762 const PetscInt s = tab->s,*is_slow,*is_fast; /* is_slow: index of slow components; is_fast: index of fast components */ 763 const PetscReal *A = tab->A,*c = tab->c; 764 PetscScalar *w = rk->work; 765 PetscScalar *y_ptr,*faststage_ptr,*slowstage_ptr; /* location to put the pointer to Y, stage_fast and stage_slow respectively */ 766 PetscScalar *ydotrhsprev_ptr,*ydotrhs_ptr; /* location to put the pointer to YdotRHS_prev and YdotRHS respectively */ 767 PetscInt i,j,k,len_slow,len_fast; /* len_slow: the number of slow comonents; len_fast: the number of fast components */ 768 PetscReal next_time_step = ts->time_step,t = ts->ptime,h = ts->time_step; 769 PetscErrorCode ierr; 770 771 PetscFunctionBegin; 772 rk->status = TS_STEP_INCOMPLETE; 773 ierr = VecDuplicate(ts->vec_sol,&stage_fast);CHKERRQ(ierr); 774 ierr = VecDuplicate(ts->vec_sol,&stage_slow);CHKERRQ(ierr); 775 ierr = VecDuplicate(ts->vec_sol,&YdotRHS_prev);CHKERRQ(ierr); 776 ierr = VecDuplicate(ts->vec_sol,&prev_sol);CHKERRQ(ierr); 777 ierr = ISGetSize(ts->iss,&len_slow);CHKERRQ(ierr); 778 ierr = ISGetSize(ts->isf,&len_fast);CHKERRQ(ierr); 779 ierr = ISGetIndices(ts->iss,&is_slow);CHKERRQ(ierr); 780 ierr = ISGetIndices(ts->isf,&is_fast);CHKERRQ(ierr); 781 ierr = ISRestoreIndices(ts->isf,&is_fast);CHKERRQ(ierr); 782 ierr = ISRestoreIndices(ts->iss,&is_slow);CHKERRQ(ierr); 783 for (i=0; i<s; i++) { 784 rk->stage_time = t + h*c[i]; 785 ierr = TSPreStage(ts,rk->stage_time);CHKERRQ(ierr); 786 ierr = VecCopy(ts->vec_sol,Y[i]);CHKERRQ(ierr); 787 for (j=0; j<i; j++) w[j] = h*A[i*s+j]; 788 ierr = VecMAXPY(Y[i],i,w,YdotRHS);CHKERRQ(ierr); 789 ierr = TSPostStage(ts,rk->stage_time,i,Y); CHKERRQ(ierr); 790 /* compute the stage RHS */ 791 ierr = TSComputeRHSFunction(ts,t+h*c[i],Y[i],YdotRHS[i]);CHKERRQ(ierr); 792 } 793 ierr = VecCopy(ts->vec_sol,prev_sol);CHKERRQ(ierr); 794 rk->slow = 0; 795 ierr = TSEvaluateStep_RKMULTIRATE(ts,tab->order,ts->vec_sol,NULL);CHKERRQ(ierr); 796 for (k=0; k<rk->dtratio; k++){ 797 for (i=0; i<s; i++) { 798 rk->stage_time = t + h/rk->dtratio*c[i]; 799 ierr = TSPreStage(ts,rk->stage_time);CHKERRQ(ierr); 800 ierr = VecCopy(prev_sol,Y[i]);CHKERRQ(ierr); 801 ierr = VecCopy(prev_sol,stage_slow);CHKERRQ(ierr); 802 ierr = VecCopy(prev_sol,stage_fast);CHKERRQ(ierr); 803 /* guarantee the stage RHS for slow components do not change while updating the stage RHS for fast components */ 804 ierr = VecCopy(YdotRHS[i],YdotRHS_prev);CHKERRQ(ierr); 805 /* update the fast components stage value */ 806 for (j=0; j<i; j++) w[j] = h/rk->dtratio*A[i*s+j]; 807 ierr = VecMAXPY(stage_fast,i,w,YdotRHS);CHKERRQ(ierr); 808 /* update the slow components stage value */ 809 ierr = VecCopy(prev_sol,Y[0]);CHKERRQ(ierr); 810 ierr = TSInterpolate_RK(ts,t+k*h/rk->dtratio+h/rk->dtratio*c[i],stage_slow);CHKERRQ(ierr); 811 /* combine the update fast components stage value and slow components stage value to Y[i] */ 812 ierr = VecGetArray(Y[i],&y_ptr);CHKERRQ(ierr); 813 ierr = VecGetArray(stage_slow,&slowstage_ptr);CHKERRQ(ierr); 814 ierr = VecGetArray(stage_fast,&faststage_ptr);CHKERRQ(ierr); 815 for (j=0; j<len_slow; j++) y_ptr[is_slow[j]] = slowstage_ptr[is_slow[j]]; 816 for (j=0; j<len_fast; j++) y_ptr[is_fast[j]] = faststage_ptr[is_fast[j]]; 817 ierr = VecRestoreArray(Y[i],&y_ptr);CHKERRQ(ierr); 818 ierr = VecRestoreArray(stage_slow,&slowstage_ptr);CHKERRQ(ierr); 819 ierr = VecRestoreArray(stage_fast,&faststage_ptr);CHKERRQ(ierr); 820 ierr = TSPostStage(ts,rk->stage_time,i,Y); CHKERRQ(ierr); 821 /* compute the stage RHS */ 822 ierr = TSComputeRHSFunction(ts,t+k*h/rk->dtratio+h/rk->dtratio*c[i],Y[i],YdotRHS[i]);CHKERRQ(ierr); 823 /* re-assign the value of slow components in YdotRHS[i] to the calculated value before working on smaller time step-size */ 824 ierr = VecGetArray(YdotRHS[i],&ydotrhs_ptr);CHKERRQ(ierr); 825 ierr = VecGetArray(YdotRHS_prev,&ydotrhsprev_ptr);CHKERRQ(ierr); 826 for (j=0; j<len_slow; j++) ydotrhs_ptr[is_slow[j]] = ydotrhsprev_ptr[is_slow[j]]; 827 ierr = VecRestoreArray(YdotRHS[i],&ydotrhs_ptr);CHKERRQ(ierr); 828 ierr = VecRestoreArray(YdotRHS_prev,&ydotrhsprev_ptr);CHKERRQ(ierr); 829 } 830 rk->slow = 1; 831 ierr = TSEvaluateStep_RKMULTIRATE(ts,tab->order,ts->vec_sol,NULL);CHKERRQ(ierr); 832 /* update the intial value for fast components */ 833 ierr = VecGetSubVector(prev_sol,ts->isf,&presol_fast);CHKERRQ(ierr); 834 ierr = VecGetSubVector(ts->vec_sol,ts->isf,&newslo_fast);CHKERRQ(ierr); 835 ierr = VecCopy(newslo_fast,presol_fast);CHKERRQ(ierr); 836 ierr = VecRestoreSubVector(prev_sol,ts->isf,&presol_fast);CHKERRQ(ierr); 837 ierr = VecRestoreSubVector(ts->vec_sol,ts->isf,&newslo_fast);CHKERRQ(ierr); 838 } 839 ts->ptime += ts->time_step; 840 ts->time_step = next_time_step; 841 rk->slow = 0; 842 rk->status = TS_STEP_COMPLETE; 843 /* free memory of work vectors */ 844 ierr = VecDestroy(&stage_fast);CHKERRQ(ierr); 845 ierr = VecDestroy(&stage_slow);CHKERRQ(ierr); 846 ierr = VecDestroy(&YdotRHS_prev);CHKERRQ(ierr); 847 ierr = VecDestroy(&prev_sol);CHKERRQ(ierr); 848 PetscFunctionReturn(0); 849 } 850 851 /* 852 This is for partitioned RHS multirate RK method 853 The step completion formula is 854 855 x1 = x0 + h b^T YdotRHS 856 857 */ 858 static PetscErrorCode TSEvaluateStep_RKPMULTIRATE(TS ts,PetscInt order,Vec X,PetscBool *done) 859 { 860 TS_RK *rk = (TS_RK*)ts->data; 861 RKTableau tab = rk->tableau; 862 Vec Xslow,Xfast; /* subvectors of X which store slow components and fast components respectively */ 863 PetscScalar *w = rk->work; 864 PetscReal h = ts->time_step; 865 PetscInt s = tab->s,j; 866 PetscErrorCode ierr; 867 868 PetscFunctionBegin; 869 ierr = VecCopy(ts->vec_sol,X);CHKERRQ(ierr); 870 if (!rk->slow){ 871 for (j=0; j<s; j++) w[j] = h*tab->b[j]; 872 ierr = VecGetSubVector(ts->vec_sol,ts->iss,&Xslow);CHKERRQ(ierr); 873 ierr = VecMAXPY(Xslow,s,w,rk->YdotRHS_slow);CHKERRQ(ierr); 874 ierr = VecRestoreSubVector(ts->vec_sol,ts->iss,&Xslow);CHKERRQ(ierr);; 875 } else {; 876 for (j=0; j<s; j++) w[j] = h/rk->dtratio*tab->b[j]; 877 ierr = VecGetSubVector(X,ts->isf,&Xfast);CHKERRQ(ierr); 878 ierr = VecMAXPY(Xfast,s,w,rk->YdotRHS_fast);CHKERRQ(ierr); 879 ierr = VecRestoreSubVector(X,ts->isf,&Xfast);CHKERRQ(ierr); 880 } 881 PetscFunctionReturn(0); 882 } 883 884 static PetscErrorCode TSStep_RKPMULTIRATE(TS ts) 885 { 886 TS_RK *rk = (TS_RK*)ts->data; 887 RKTableau tab = rk->tableau; 888 Vec *Y = rk->Y; 889 Vec *YdotRHS_fast = rk->YdotRHS_fast,*YdotRHS_slow = rk->YdotRHS_slow; 890 Vec Yslow,Yfast; /* subvectors store the stges of slow components and fast components respectively */ 891 Vec Yfast_prev,Yfast_new,prev_sol; /* subvectors store the previous and new solution of fast components and vector store the previous solution */ 892 const PetscInt s = tab->s; 893 const PetscReal *A = tab->A,*c = tab->c; 894 PetscScalar *w = rk->work; 895 PetscInt i,j,k; 896 PetscReal next_time_step = ts->time_step,t = ts->ptime,h = ts->time_step; 897 PetscErrorCode ierr; 898 899 PetscFunctionBegin; 900 rk->status = TS_STEP_INCOMPLETE; 901 for (i=0; i<s; i++) { 902 rk->stage_time = t + h*c[i]; 903 ierr = TSPreStage(ts,rk->stage_time);CHKERRQ(ierr); 904 ierr = VecCopy(ts->vec_sol,Y[i]);CHKERRQ(ierr); 905 /*ierr = VecCopy(ts->vec_sol,Yc[i]);CHKERRQ(ierr);*/ 906 ierr = VecGetSubVector(Y[i],ts->iss,&Yslow);CHKERRQ(ierr); 907 ierr = VecGetSubVector(Y[i],ts->isf,&Yfast);CHKERRQ(ierr); 908 for (j=0; j<i; j++) w[j] = h*A[i*s+j]; 909 ierr = VecMAXPY(Yslow,i,w,YdotRHS_slow);CHKERRQ(ierr); 910 ierr = VecMAXPY(Yfast,i,w,YdotRHS_fast);CHKERRQ(ierr); 911 ierr = VecRestoreSubVector(Y[i],ts->iss,&Yslow);CHKERRQ(ierr); 912 ierr = VecRestoreSubVector(Y[i],ts->isf,&Yfast);CHKERRQ(ierr); 913 ierr = TSPostStage(ts,rk->stage_time,i,Y); CHKERRQ(ierr); 914 /* compute the stage RHS for both slow and fast components */ 915 ierr = TSComputeRHSFunctionslow(ts,t+h*c[i],Y[i],YdotRHS_slow[i]);CHKERRQ(ierr); 916 ierr = TSComputeRHSFunctionfast(ts,t+h*c[i],Y[i],YdotRHS_fast[i]);CHKERRQ(ierr); 917 } 918 /* store the value of slow components at previous time */ 919 ierr = VecDuplicate(ts->vec_sol,&prev_sol);CHKERRQ(ierr); 920 ierr = VecCopy(ts->vec_sol,prev_sol);CHKERRQ(ierr); 921 rk->slow = 0; 922 ierr = TSEvaluateStep_RKPMULTIRATE(ts,tab->order,ts->vec_sol,NULL);CHKERRQ(ierr); 923 for (k=0; k<rk->dtratio; k++){ 924 for (i=0; i<s; i++) { 925 rk->stage_time = t + h/rk->dtratio*c[i]; 926 ierr = TSPreStage(ts,rk->stage_time);CHKERRQ(ierr); 927 ierr = VecCopy(prev_sol,Y[i]);CHKERRQ(ierr); 928 /* stage value for fast components */ 929 for (j=0; j<i; j++) w[j] = h/rk->dtratio*A[i*s+j]; 930 ierr = VecGetSubVector(Y[i],ts->isf,&Yfast);CHKERRQ(ierr); 931 ierr = VecMAXPY(Yfast,i,w,YdotRHS_fast);CHKERRQ(ierr); 932 ierr = VecRestoreSubVector(Y[i],ts->isf,&Yfast);CHKERRQ(ierr); 933 /* stage value for slow components */ 934 ierr = VecCopy(prev_sol,Y[0]);CHKERRQ(ierr); 935 ierr = VecGetSubVector(Y[i],ts->iss,&Yslow);CHKERRQ(ierr); 936 ierr = TSInterpolate_PARTITIONEDMRK(ts,t+k*h/rk->dtratio+h/rk->dtratio*c[i],Yslow);CHKERRQ(ierr); 937 ierr = VecRestoreSubVector(Y[i],ts->iss,&Yslow);CHKERRQ(ierr); 938 ierr = TSPostStage(ts,rk->stage_time,i,Y); CHKERRQ(ierr); 939 /* compute the stage RHS for fast components */ 940 ierr = TSComputeRHSFunctionfast(ts,t+k*h/rk->dtratio+h/rk->dtratio*c[i],Y[i],YdotRHS_fast[i]);CHKERRQ(ierr); 941 } 942 /* update the value of fast components whenusing fast time step */ 943 rk->slow = 1; 944 ierr = TSEvaluateStep_RKPMULTIRATE(ts,tab->order,ts->vec_sol,NULL);CHKERRQ(ierr); 945 ierr = VecGetSubVector(prev_sol,ts->isf,&Yfast_prev);CHKERRQ(ierr); 946 ierr = VecGetSubVector(ts->vec_sol,ts->isf,&Yfast_new);CHKERRQ(ierr); 947 ierr = VecCopy(Yfast_new,Yfast_prev);CHKERRQ(ierr); 948 ierr = VecRestoreSubVector(prev_sol,ts->isf,&Yfast_prev);CHKERRQ(ierr); 949 ierr = VecRestoreSubVector(ts->vec_sol,ts->isf,&Yfast_new);CHKERRQ(ierr); 950 } 951 ts->ptime += ts->time_step; 952 ts->time_step = next_time_step; 953 rk->slow = 0; 954 rk->status = TS_STEP_COMPLETE; 955 ierr = VecDestroy(&prev_sol);CHKERRQ(ierr); 956 PetscFunctionReturn(0); 957 } 958 959 static PetscErrorCode TSAdjointSetUp_RK(TS ts) 960 { 961 TS_RK *rk = (TS_RK*)ts->data; 962 RKTableau tab = rk->tableau; 963 PetscInt s = tab->s; 964 PetscErrorCode ierr; 965 966 PetscFunctionBegin; 967 if (ts->adjointsetupcalled++) PetscFunctionReturn(0); 968 ierr = VecDuplicateVecs(ts->vecs_sensi[0],s*ts->numcost,&rk->VecDeltaLam);CHKERRQ(ierr); 969 if(ts->vecs_sensip) { 970 ierr = VecDuplicateVecs(ts->vecs_sensip[0],s*ts->numcost,&rk->VecDeltaMu);CHKERRQ(ierr); 971 } 972 PetscFunctionReturn(0); 973 } 974 975 static PetscErrorCode TSAdjointStep_RK(TS ts) 976 { 977 TS_RK *rk = (TS_RK*)ts->data; 978 RKTableau tab = rk->tableau; 979 const PetscInt s = tab->s; 980 const PetscReal *A = tab->A,*b = tab->b,*c = tab->c; 981 PetscScalar *w = rk->work; 982 Vec *Y = rk->Y,*VecDeltaLam = rk->VecDeltaLam,*VecDeltaMu = rk->VecDeltaMu; 983 PetscInt i,j,nadj; 984 PetscReal t = ts->ptime; 985 PetscErrorCode ierr; 986 PetscReal h = ts->time_step; 987 988 PetscFunctionBegin; 989 rk->status = TS_STEP_INCOMPLETE; 990 for (i=s-1; i>=0; i--) { 991 Mat J; 992 PetscReal stage_time = t + h*(1.0-c[i]); 993 PetscBool zero = PETSC_FALSE; 994 995 ierr = TSGetRHSJacobian(ts,&J,NULL,NULL,NULL);CHKERRQ(ierr); 996 ierr = TSComputeRHSJacobian(ts,stage_time,Y[i],J,J);CHKERRQ(ierr); 997 if (ts->vec_costintegral) { 998 ierr = TSAdjointComputeDRDYFunction(ts,stage_time,Y[i],ts->vecs_drdy);CHKERRQ(ierr); 999 } 1000 /* Stage values of mu */ 1001 if (ts->vecs_sensip) { 1002 ierr = TSAdjointComputeRHSJacobian(ts,stage_time,Y[i],ts->Jacp);CHKERRQ(ierr); 1003 if (ts->vec_costintegral) { 1004 ierr = TSAdjointComputeDRDPFunction(ts,stage_time,Y[i],ts->vecs_drdp);CHKERRQ(ierr); 1005 } 1006 } 1007 1008 if (b[i] == 0 && i == s-1) zero = PETSC_TRUE; 1009 for (nadj=0; nadj<ts->numcost; nadj++) { 1010 DM dm; 1011 Vec VecSensiTemp; 1012 1013 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 1014 ierr = DMGetGlobalVector(dm,&VecSensiTemp);CHKERRQ(ierr); 1015 /* Stage values of lambda */ 1016 if (!zero) { 1017 ierr = VecCopy(ts->vecs_sensi[nadj],VecSensiTemp);CHKERRQ(ierr); 1018 ierr = VecScale(VecSensiTemp,-h*b[i]);CHKERRQ(ierr); 1019 for (j=i+1; j<s; j++) w[j-i-1] = -h*A[j*s+i]; 1020 ierr = VecMAXPY(VecSensiTemp,s-i-1,w,&VecDeltaLam[nadj*s+i+1]);CHKERRQ(ierr); 1021 ierr = MatMultTranspose(J,VecSensiTemp,VecDeltaLam[nadj*s+i]);CHKERRQ(ierr); 1022 } else { 1023 ierr = VecSet(VecDeltaLam[nadj*s+i],0);CHKERRQ(ierr); 1024 } 1025 if (ts->vec_costintegral) { 1026 ierr = VecAXPY(VecDeltaLam[nadj*s+i],-h*b[i],ts->vecs_drdy[nadj]);CHKERRQ(ierr); 1027 } 1028 1029 /* Stage values of mu */ 1030 if (ts->vecs_sensip) { 1031 if (!zero) { 1032 ierr = MatMultTranspose(ts->Jacp,VecSensiTemp,VecDeltaMu[nadj*s+i]);CHKERRQ(ierr); 1033 } else { 1034 ierr = VecSet(VecDeltaMu[nadj*s+i],0);CHKERRQ(ierr); 1035 } 1036 if (ts->vec_costintegral) { 1037 ierr = VecAXPY(VecDeltaMu[nadj*s+i],-h*b[i],ts->vecs_drdp[nadj]);CHKERRQ(ierr); 1038 } 1039 } 1040 ierr = DMRestoreGlobalVector(dm,&VecSensiTemp);CHKERRQ(ierr); 1041 } 1042 } 1043 1044 for (j=0; j<s; j++) w[j] = 1.0; 1045 for (nadj=0; nadj<ts->numcost; nadj++) { 1046 ierr = VecMAXPY(ts->vecs_sensi[nadj],s,w,&VecDeltaLam[nadj*s]);CHKERRQ(ierr); 1047 if (ts->vecs_sensip) { 1048 ierr = VecMAXPY(ts->vecs_sensip[nadj],s,w,&VecDeltaMu[nadj*s]);CHKERRQ(ierr); 1049 } 1050 } 1051 rk->status = TS_STEP_COMPLETE; 1052 PetscFunctionReturn(0); 1053 } 1054 1055 static PetscErrorCode TSRKTableauReset(TS ts) 1056 { 1057 TS_RK *rk = (TS_RK*)ts->data; 1058 RKTableau tab = rk->tableau; 1059 PetscErrorCode ierr; 1060 1061 PetscFunctionBegin; 1062 if (!tab) PetscFunctionReturn(0); 1063 ierr = PetscFree(rk->work);CHKERRQ(ierr); 1064 ierr = VecDestroyVecs(tab->s,&rk->Y);CHKERRQ(ierr); 1065 ierr = VecDestroyVecs(tab->s,&rk->YdotRHS);CHKERRQ(ierr); 1066 ierr = VecDestroyVecs(tab->s,&rk->YdotRHS_fast);CHKERRQ(ierr); 1067 ierr = VecDestroyVecs(tab->s,&rk->YdotRHS_slow);CHKERRQ(ierr); 1068 ierr = VecDestroyVecs(tab->s*ts->numcost,&rk->VecDeltaLam);CHKERRQ(ierr); 1069 ierr = VecDestroyVecs(tab->s*ts->numcost,&rk->VecDeltaMu);CHKERRQ(ierr); 1070 PetscFunctionReturn(0); 1071 } 1072 1073 static PetscErrorCode TSReset_RK(TS ts) 1074 { 1075 TS_RK *rk = (TS_RK*)ts->data; 1076 PetscErrorCode ierr; 1077 1078 PetscFunctionBegin; 1079 ierr = TSRKTableauReset(ts);CHKERRQ(ierr); 1080 ierr = VecDestroy(&rk->VecCostIntegral0);CHKERRQ(ierr); 1081 PetscFunctionReturn(0); 1082 } 1083 1084 static PetscErrorCode DMCoarsenHook_TSRK(DM fine,DM coarse,void *ctx) 1085 { 1086 PetscFunctionBegin; 1087 PetscFunctionReturn(0); 1088 } 1089 1090 static PetscErrorCode DMRestrictHook_TSRK(DM fine,Mat restrct,Vec rscale,Mat inject,DM coarse,void *ctx) 1091 { 1092 PetscFunctionBegin; 1093 PetscFunctionReturn(0); 1094 } 1095 1096 1097 static PetscErrorCode DMSubDomainHook_TSRK(DM dm,DM subdm,void *ctx) 1098 { 1099 PetscFunctionBegin; 1100 PetscFunctionReturn(0); 1101 } 1102 1103 static PetscErrorCode DMSubDomainRestrictHook_TSRK(DM dm,VecScatter gscat,VecScatter lscat,DM subdm,void *ctx) 1104 { 1105 1106 PetscFunctionBegin; 1107 PetscFunctionReturn(0); 1108 } 1109 /* 1110 static PetscErrorCode RKSetAdjCoe(RKTableau tab) 1111 { 1112 PetscReal *A,*b; 1113 PetscInt s,i,j; 1114 PetscErrorCode ierr; 1115 1116 PetscFunctionBegin; 1117 s = tab->s; 1118 ierr = PetscMalloc2(s*s,&A,s,&b);CHKERRQ(ierr); 1119 1120 for (i=0; i<s; i++) 1121 for (j=0; j<s; j++) { 1122 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]; 1123 ierr = PetscPrintf(PETSC_COMM_WORLD,"Coefficients: A[%D][%D]=%.6f\n",i,j,A[i*s+j]);CHKERRQ(ierr); 1124 } 1125 1126 for (i=0; i<s; i++) b[i] = (tab->b[s-1-i]==0)? 0: -tab->b[s-1-i]; 1127 1128 ierr = PetscMemcpy(tab->A,A,s*s*sizeof(A[0]));CHKERRQ(ierr); 1129 ierr = PetscMemcpy(tab->b,b,s*sizeof(b[0]));CHKERRQ(ierr); 1130 ierr = PetscFree2(A,b);CHKERRQ(ierr); 1131 PetscFunctionReturn(0); 1132 } 1133 */ 1134 1135 static PetscErrorCode TSRKTableauSetUp(TS ts) 1136 { 1137 TS_RK *rk = (TS_RK*)ts->data; 1138 RKTableau tab = rk->tableau; 1139 Vec YdotRHS_fast,YdotRHS_slow; 1140 PetscErrorCode ierr; 1141 1142 PetscFunctionBegin; 1143 ierr = PetscMalloc1(tab->s,&rk->work);CHKERRQ(ierr); 1144 ierr = VecDuplicateVecs(ts->vec_sol,tab->s,&rk->Y);CHKERRQ(ierr); 1145 ierr = VecDuplicateVecs(ts->vec_sol,tab->s,&rk->YdotRHS);CHKERRQ(ierr); 1146 ierr = VecGetSubVector(ts->vec_sol,ts->isf,&YdotRHS_fast);CHKERRQ(ierr); 1147 ierr = VecDuplicateVecs(YdotRHS_fast,tab->s,&rk->YdotRHS_fast);CHKERRQ(ierr); 1148 ierr = VecRestoreSubVector(ts->vec_sol,ts->isf,&YdotRHS_fast);CHKERRQ(ierr); 1149 ierr = VecGetSubVector(ts->vec_sol,ts->iss,&YdotRHS_slow);CHKERRQ(ierr); 1150 ierr = VecDuplicateVecs(YdotRHS_slow,tab->s,&rk->YdotRHS_slow);CHKERRQ(ierr); 1151 ierr = VecRestoreSubVector(ts->vec_sol,ts->iss,&YdotRHS_slow);CHKERRQ(ierr); 1152 1153 PetscFunctionReturn(0); 1154 } 1155 1156 1157 static PetscErrorCode TSSetUp_RK(TS ts) 1158 { 1159 TS_RK *rk = (TS_RK*)ts->data; 1160 PetscErrorCode ierr; 1161 DM dm; 1162 1163 PetscFunctionBegin; 1164 ierr = TSCheckImplicitTerm(ts);CHKERRQ(ierr); 1165 ierr = TSRKTableauSetUp(ts);CHKERRQ(ierr); 1166 if (!rk->VecCostIntegral0 && ts->vec_costintegral && ts->costintegralfwd) { /* back up cost integral */ 1167 ierr = VecDuplicate(ts->vec_costintegral,&rk->VecCostIntegral0);CHKERRQ(ierr); 1168 } 1169 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 1170 ierr = DMCoarsenHookAdd(dm,DMCoarsenHook_TSRK,DMRestrictHook_TSRK,ts);CHKERRQ(ierr); 1171 ierr = DMSubDomainHookAdd(dm,DMSubDomainHook_TSRK,DMSubDomainRestrictHook_TSRK,ts);CHKERRQ(ierr); 1172 PetscFunctionReturn(0); 1173 } 1174 1175 /* 1176 construct a database to choose single-step rk method, or combined rhs mutirate rk method or partitioned rhs rk method 1177 */ 1178 const char *const TSRKMultirateTypes[] = {"NONE","COMBINED","PARTITIONED","TSRKMultirateType","RKM_",0}; 1179 1180 static PetscErrorCode TSSetFromOptions_RK(PetscOptionItems *PetscOptionsObject,TS ts) 1181 { 1182 TS_RK *rk = (TS_RK*)ts->data; 1183 PetscErrorCode ierr; 1184 1185 PetscFunctionBegin; 1186 ierr = PetscOptionsHead(PetscOptionsObject,"RK ODE solver options");CHKERRQ(ierr); 1187 { 1188 RKTableauLink link; 1189 PetscInt count,choice; 1190 PetscBool flg; 1191 const char **namelist; 1192 PetscInt rkmtype = 0; 1193 1194 for (link=RKTableauList,count=0; link; link=link->next,count++) ; 1195 ierr = PetscMalloc1(count,(char***)&namelist);CHKERRQ(ierr); 1196 for (link=RKTableauList,count=0; link; link=link->next,count++) namelist[count] = link->tab.name; 1197 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); 1198 if (flg) {ierr = TSRKSetMultirateType(ts,rkmtype);CHKERRQ(ierr);} 1199 ierr = PetscOptionsEList("-ts_rk_type","Family of RK method","TSRKSetType",(const char*const*)namelist,count,rk->tableau->name,&choice,&flg);CHKERRQ(ierr); 1200 if (flg) {ierr = TSRKSetType(ts,namelist[choice]);CHKERRQ(ierr);} 1201 ierr = PetscFree(namelist);CHKERRQ(ierr); 1202 } 1203 ierr = PetscOptionsTail();CHKERRQ(ierr); 1204 ierr = PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"Multirate methods options","");CHKERRQ(ierr); 1205 ierr = PetscOptionsInt("-ts_rk_dtratio","time step ratio between slow and fast","",rk->dtratio,&rk->dtratio,NULL);CHKERRQ(ierr); 1206 ierr = PetscOptionsEnd();CHKERRQ(ierr); 1207 PetscFunctionReturn(0); 1208 } 1209 1210 static PetscErrorCode TSView_RK(TS ts,PetscViewer viewer) 1211 { 1212 TS_RK *rk = (TS_RK*)ts->data; 1213 PetscBool iascii; 1214 PetscErrorCode ierr; 1215 1216 PetscFunctionBegin; 1217 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 1218 if (iascii) { 1219 RKTableau tab = rk->tableau; 1220 TSRKType rktype; 1221 char buf[512]; 1222 ierr = TSRKGetType(ts,&rktype);CHKERRQ(ierr); 1223 ierr = PetscViewerASCIIPrintf(viewer," RK type %s\n",rktype);CHKERRQ(ierr); 1224 ierr = PetscViewerASCIIPrintf(viewer," Order: %D\n",tab->order);CHKERRQ(ierr); 1225 ierr = PetscViewerASCIIPrintf(viewer," FSAL property: %s\n",tab->FSAL ? "yes" : "no");CHKERRQ(ierr); 1226 ierr = PetscFormatRealArray(buf,sizeof(buf),"% 8.6f",tab->s,tab->c);CHKERRQ(ierr); 1227 ierr = PetscViewerASCIIPrintf(viewer," Abscissa c = %s\n",buf);CHKERRQ(ierr); 1228 } 1229 PetscFunctionReturn(0); 1230 } 1231 1232 static PetscErrorCode TSLoad_RK(TS ts,PetscViewer viewer) 1233 { 1234 PetscErrorCode ierr; 1235 TSAdapt adapt; 1236 1237 PetscFunctionBegin; 1238 ierr = TSGetAdapt(ts,&adapt);CHKERRQ(ierr); 1239 ierr = TSAdaptLoad(adapt,viewer);CHKERRQ(ierr); 1240 PetscFunctionReturn(0); 1241 } 1242 1243 /*@C 1244 TSRKSetType - Set the type of RK scheme 1245 1246 Logically collective 1247 1248 Input Parameter: 1249 + ts - timestepping context 1250 - rktype - type of RK-scheme 1251 1252 Options Database: 1253 . -ts_rk_type - <1fe,2a,3,3bs,4,5f,5dp,5bs> 1254 1255 Level: intermediate 1256 1257 .seealso: TSRKGetType(), TSRK, TSRKType, TSRK1FE, TSRK2A, TSRK3, TSRK3BS, TSRK4, TSRK5F, TSRK5DP, TSRK5BS 1258 @*/ 1259 PetscErrorCode TSRKSetType(TS ts,TSRKType rktype) 1260 { 1261 PetscErrorCode ierr; 1262 1263 PetscFunctionBegin; 1264 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1265 PetscValidCharPointer(rktype,2); 1266 ierr = PetscTryMethod(ts,"TSRKSetType_C",(TS,TSRKType),(ts,rktype));CHKERRQ(ierr); 1267 PetscFunctionReturn(0); 1268 } 1269 1270 /*@C 1271 TSRKSetMultirateType - Set the type of RK Multirate scheme 1272 1273 Logically collective 1274 1275 Input Parameter: 1276 + ts - timestepping context 1277 - rkmtype - type of RKM-scheme 1278 1279 Options Database: 1280 . -ts_rk_multiarte_type - <none,combined,partitioned> 1281 1282 Level: intermediate 1283 @*/ 1284 PetscErrorCode TSRKSetMultirateType(TS ts, TSRKMultirateType rkmtype) 1285 { 1286 TS_RK *rk = (TS_RK*)ts->data; 1287 PetscErrorCode ierr; 1288 1289 PetscFunctionBegin; 1290 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1291 switch(rkmtype){ 1292 case RKM_NONE: 1293 break; 1294 case RKM_COMBINED: 1295 ts->ops->step = TSStep_RKMULTIRATE; 1296 ts->ops->evaluatestep = TSEvaluateStep_RKMULTIRATE; 1297 rk->dtratio = 2; 1298 ierr = TSRKSetType(ts,TSMRKDefault);CHKERRQ(ierr); 1299 break; 1300 case RKM_PARTITIONED: 1301 ts->ops->step = TSStep_RKPMULTIRATE; 1302 ts->ops->evaluatestep = TSEvaluateStep_RKPMULTIRATE; 1303 ts->ops->interpolate = TSInterpolate_PARTITIONEDMRK; 1304 rk->dtratio = 2; 1305 ierr = TSRKSetType(ts,TSMRKDefault);CHKERRQ(ierr); 1306 break; 1307 default : 1308 SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown type '%s'",rkmtype); 1309 } 1310 PetscFunctionReturn(0); 1311 } 1312 1313 /*@C 1314 TSRKGetType - Get the type of RK scheme 1315 1316 Logically collective 1317 1318 Input Parameter: 1319 . ts - timestepping context 1320 1321 Output Parameter: 1322 . rktype - type of RK-scheme 1323 1324 Level: intermediate 1325 1326 .seealso: TSRKGetType() 1327 @*/ 1328 PetscErrorCode TSRKGetType(TS ts,TSRKType *rktype) 1329 { 1330 PetscErrorCode ierr; 1331 1332 PetscFunctionBegin; 1333 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1334 ierr = PetscUseMethod(ts,"TSRKGetType_C",(TS,TSRKType*),(ts,rktype));CHKERRQ(ierr); 1335 PetscFunctionReturn(0); 1336 } 1337 1338 static PetscErrorCode TSRKGetType_RK(TS ts,TSRKType *rktype) 1339 { 1340 TS_RK *rk = (TS_RK*)ts->data; 1341 1342 PetscFunctionBegin; 1343 *rktype = rk->tableau->name; 1344 PetscFunctionReturn(0); 1345 } 1346 1347 static PetscErrorCode TSRKSetType_RK(TS ts,TSRKType rktype) 1348 { 1349 TS_RK *rk = (TS_RK*)ts->data; 1350 PetscErrorCode ierr; 1351 PetscBool match; 1352 RKTableauLink link; 1353 1354 PetscFunctionBegin; 1355 if (rk->tableau) { 1356 ierr = PetscStrcmp(rk->tableau->name,rktype,&match);CHKERRQ(ierr); 1357 if (match) PetscFunctionReturn(0); 1358 } 1359 for (link = RKTableauList; link; link=link->next) { 1360 ierr = PetscStrcmp(link->tab.name,rktype,&match);CHKERRQ(ierr); 1361 if (match) { 1362 if (ts->setupcalled) {ierr = TSRKTableauReset(ts);CHKERRQ(ierr);} 1363 rk->tableau = &link->tab; 1364 if (ts->setupcalled) {ierr = TSRKTableauSetUp(ts);CHKERRQ(ierr);} 1365 ts->default_adapt_type = rk->tableau->bembed ? TSADAPTBASIC : TSADAPTNONE; 1366 PetscFunctionReturn(0); 1367 } 1368 } 1369 SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_UNKNOWN_TYPE,"Could not find '%s'",rktype); 1370 PetscFunctionReturn(0); 1371 } 1372 1373 static PetscErrorCode TSGetStages_RK(TS ts,PetscInt *ns,Vec **Y) 1374 { 1375 TS_RK *rk = (TS_RK*)ts->data; 1376 1377 PetscFunctionBegin; 1378 if (ns) *ns = rk->tableau->s; 1379 if (Y) *Y = rk->Y; 1380 PetscFunctionReturn(0); 1381 } 1382 1383 static PetscErrorCode TSDestroy_RK(TS ts) 1384 { 1385 PetscErrorCode ierr; 1386 1387 PetscFunctionBegin; 1388 ierr = TSReset_RK(ts);CHKERRQ(ierr); 1389 if (ts->dm) { 1390 ierr = DMCoarsenHookRemove(ts->dm,DMCoarsenHook_TSRK,DMRestrictHook_TSRK,ts);CHKERRQ(ierr); 1391 ierr = DMSubDomainHookRemove(ts->dm,DMSubDomainHook_TSRK,DMSubDomainRestrictHook_TSRK,ts);CHKERRQ(ierr); 1392 } 1393 ierr = PetscFree(ts->data);CHKERRQ(ierr); 1394 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSRKGetType_C",NULL);CHKERRQ(ierr); 1395 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSRKSetType_C",NULL);CHKERRQ(ierr); 1396 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSRKSetMultirateType_C",NULL);CHKERRQ(ierr); 1397 PetscFunctionReturn(0); 1398 } 1399 1400 /*MC 1401 TSRK - ODE and DAE solver using Runge-Kutta schemes 1402 1403 The user should provide the right hand side of the equation 1404 using TSSetRHSFunction(). 1405 1406 Notes: 1407 The default is TSRK3BS, it can be changed with TSRKSetType() or -ts_rk_type 1408 1409 Level: beginner 1410 1411 .seealso: TSCreate(), TS, TSSetType(), TSRKSetType(), TSRKGetType(), TSRKSetFullyImplicit(), TSRK2D, TTSRK2E, TSRK3, 1412 TSRK4, TSRK5, TSRKPRSSP2, TSRKBPR3, TSRKType, TSRKRegister(), TSRKSetMultirateType() 1413 1414 M*/ 1415 PETSC_EXTERN PetscErrorCode TSCreate_RK(TS ts) 1416 { 1417 TS_RK *rk; 1418 PetscErrorCode ierr; 1419 1420 PetscFunctionBegin; 1421 ierr = TSRKInitializePackage();CHKERRQ(ierr); 1422 1423 ts->ops->reset = TSReset_RK; 1424 ts->ops->destroy = TSDestroy_RK; 1425 ts->ops->view = TSView_RK; 1426 ts->ops->load = TSLoad_RK; 1427 ts->ops->setup = TSSetUp_RK; 1428 ts->ops->adjointsetup = TSAdjointSetUp_RK; 1429 ts->ops->interpolate = TSInterpolate_RK; 1430 ts->ops->step = TSStep_RK; 1431 ts->ops->evaluatestep = TSEvaluateStep_RK; 1432 ts->ops->rollback = TSRollBack_RK; 1433 ts->ops->setfromoptions = TSSetFromOptions_RK; 1434 ts->ops->getstages = TSGetStages_RK; 1435 ts->ops->adjointstep = TSAdjointStep_RK; 1436 1437 ts->ops->adjointintegral = TSAdjointCostIntegral_RK; 1438 ts->ops->forwardintegral = TSForwardCostIntegral_RK; 1439 1440 ierr = PetscNewLog(ts,&rk);CHKERRQ(ierr); 1441 ts->data = (void*)rk; 1442 1443 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSRKGetType_C",TSRKGetType_RK);CHKERRQ(ierr); 1444 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSRKSetType_C",TSRKSetType_RK);CHKERRQ(ierr); 1445 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSRKSetMultirateType_C",TSRKSetMultirateType);CHKERRQ(ierr); 1446 1447 ierr = TSRKSetType(ts,TSRKDefault);CHKERRQ(ierr); 1448 PetscFunctionReturn(0); 1449 } 1450