1 /* 2 Code for timestepping with Rosenbrock W methods 3 4 Notes: 5 The general system is written as 6 7 F(t,U,Udot) = G(t,U) 8 9 where F represents the stiff part of the physics and G represents the non-stiff part. 10 This method is designed to be linearly implicit on F and can use an approximate and lagged Jacobian. 11 12 */ 13 #include <petsc/private/tsimpl.h> /*I "petscts.h" I*/ 14 #include <petscdm.h> 15 16 #include <petsc/private/kernels/blockinvert.h> 17 18 static TSRosWType TSRosWDefault = TSROSWRA34PW2; 19 static PetscBool TSRosWRegisterAllCalled; 20 static PetscBool TSRosWPackageInitialized; 21 22 typedef struct _RosWTableau *RosWTableau; 23 struct _RosWTableau { 24 char *name; 25 PetscInt order; /* Classical approximation order of the method */ 26 PetscInt s; /* Number of stages */ 27 PetscInt pinterp; /* Interpolation order */ 28 PetscReal *A; /* Propagation table, strictly lower triangular */ 29 PetscReal *Gamma; /* Stage table, lower triangular with nonzero diagonal */ 30 PetscBool *GammaZeroDiag; /* Diagonal entries that are zero in stage table Gamma, vector indicating explicit statages */ 31 PetscReal *GammaExplicitCorr; /* Coefficients for correction terms needed for explicit stages in transformed variables*/ 32 PetscReal *b; /* Step completion table */ 33 PetscReal *bembed; /* Step completion table for embedded method of order one less */ 34 PetscReal *ASum; /* Row sum of A */ 35 PetscReal *GammaSum; /* Row sum of Gamma, only needed for non-autonomous systems */ 36 PetscReal *At; /* Propagation table in transformed variables */ 37 PetscReal *bt; /* Step completion table in transformed variables */ 38 PetscReal *bembedt; /* Step completion table of order one less in transformed variables */ 39 PetscReal *GammaInv; /* Inverse of Gamma, used for transformed variables */ 40 PetscReal ccfl; /* Placeholder for CFL coefficient relative to forward Euler */ 41 PetscReal *binterpt; /* Dense output formula */ 42 }; 43 typedef struct _RosWTableauLink *RosWTableauLink; 44 struct _RosWTableauLink { 45 struct _RosWTableau tab; 46 RosWTableauLink next; 47 }; 48 static RosWTableauLink RosWTableauList; 49 50 typedef struct { 51 RosWTableau tableau; 52 Vec *Y; /* States computed during the step, used to complete the step */ 53 Vec Ydot; /* Work vector holding Ydot during residual evaluation */ 54 Vec Ystage; /* Work vector for the state value at each stage */ 55 Vec Zdot; /* Ydot = Zdot + shift*Y */ 56 Vec Zstage; /* Y = Zstage + Y */ 57 Vec vec_sol_prev; /* Solution from the previous step (used for interpolation and rollback)*/ 58 PetscScalar *work; /* Scalar work space of length number of stages, used to prepare VecMAXPY() */ 59 PetscReal scoeff; /* shift = scoeff/dt */ 60 PetscReal stage_time; 61 PetscReal stage_explicit; /* Flag indicates that the current stage is explicit */ 62 PetscBool recompute_jacobian; /* Recompute the Jacobian at each stage, default is to freeze the Jacobian at the start of each step */ 63 TSStepStatus status; 64 } TS_RosW; 65 66 /*MC 67 TSROSWTHETA1 - One stage first order L-stable Rosenbrock-W scheme (aka theta method). 68 69 Only an approximate Jacobian is needed. 70 71 Level: intermediate 72 73 .seealso: [](ch_ts), `TSROSW` 74 M*/ 75 76 /*MC 77 TSROSWTHETA2 - One stage second order A-stable Rosenbrock-W scheme (aka theta method). 78 79 Only an approximate Jacobian is needed. 80 81 Level: intermediate 82 83 .seealso: [](ch_ts), `TSROSW` 84 M*/ 85 86 /*MC 87 TSROSW2M - Two stage second order L-stable Rosenbrock-W scheme. 88 89 Only an approximate Jacobian is needed. By default, it is only recomputed once per step. This method is a reflection of TSROSW2P. 90 91 Level: intermediate 92 93 .seealso: [](ch_ts), `TSROSW` 94 M*/ 95 96 /*MC 97 TSROSW2P - Two stage second order L-stable Rosenbrock-W scheme. 98 99 Only an approximate Jacobian is needed. By default, it is only recomputed once per step. This method is a reflection of TSROSW2M. 100 101 Level: intermediate 102 103 .seealso: [](ch_ts), `TSROSW` 104 M*/ 105 106 /*MC 107 TSROSWRA3PW - Three stage third order Rosenbrock-W scheme for PDAE of index 1 {cite}`rang_2005` 108 109 Only an approximate Jacobian is needed. By default, it is only recomputed once per step. 110 111 This is strongly A-stable with R(infty) = 0.73. The embedded method of order 2 is strongly A-stable with R(infty) = 0.73. 112 113 Level: intermediate 114 115 .seealso: [](ch_ts), `TSROSW` 116 M*/ 117 118 /*MC 119 TSROSWRA34PW2 - Four stage third order L-stable Rosenbrock-W scheme for PDAE of index 1 {cite}`rang_2005`. 120 121 Only an approximate Jacobian is needed. By default, it is only recomputed once per step. 122 123 This is strongly A-stable with R(infty) = 0. The embedded method of order 2 is strongly A-stable with R(infty) = 0.48. 124 125 Level: intermediate 126 127 .seealso: [](ch_ts), `TSROSW` 128 M*/ 129 130 /*MC 131 TSROSWR34PRW - Four stage third order L-stable Rosenbrock-W scheme for PDAE of index 1 {cite}`rang2015improved`. 132 133 Only an approximate Jacobian is needed. By default, it is only recomputed once per step. 134 135 This is strongly A-stable with R(infty) = 0. The embedded method of order 2 is strongly A-stable with R(infty) = 0.25. 136 This method is B_{PR} consistent of order 3. 137 This method is spelled "ROS34PRw" in the paper, an improvement to an earlier "ROS34PRW" method from the same author with B_{PR} order 2. 138 139 Level: intermediate 140 141 .seealso: [](ch_ts), `TSROSW` 142 M*/ 143 144 /*MC 145 TSROSWR3PRL2 - Four stage third order L-stable Rosenbrock-W scheme for PDAE of index 1 {cite}`rang2015improved`. 146 147 Only an approximate Jacobian is needed. By default, it is only recomputed once per step. 148 149 This is strongly A-stable with R(infty) = 0. The embedded method of order 2 is strongly A-stable with R(infty) = 0.25. 150 This method is B_{PR} consistent of order 3. 151 152 Level: intermediate 153 154 .seealso: [](ch_ts), `TSROSW` 155 M*/ 156 157 /*MC 158 TSROSWRODAS3 - Four stage third order L-stable Rosenbrock scheme {cite}`sandu_1997` 159 160 By default, the Jacobian is only recomputed once per step. 161 162 Both the third order and embedded second order methods are stiffly accurate and L-stable. 163 164 Level: intermediate 165 166 .seealso: [](ch_ts), `TSROSW`, `TSROSWSANDU3` 167 M*/ 168 169 /*MC 170 TSROSWRODASPR - Six stage fourth order L-stable Rosenbrock scheme {cite}`rang2015improved` 171 172 By default, the Jacobian is only recomputed once per step. 173 174 Both the fourth order and embedded third order methods are stiffly accurate and L-stable. 175 The method is B_{PR} consistent of order 3, which ensures convergence order for non-stiff, medium stiff, and stiff problems. 176 177 Level: intermediate 178 179 .seealso: [](ch_ts), `TSROSW`, `TSROSWR34PRW`, `TSROSWR3PRL2` 180 M*/ 181 182 /*MC 183 TSROSWRODASPR2 - Six stage fourth order L-stable Rosenbrock scheme {cite}`rang2015improved` 184 185 By default, the Jacobian is only recomputed once per step. 186 187 Both the fourth order and embedded third order methods are stiffly accurate and L-stable. 188 The method is B_{PR} consistent of order 3, which ensures convergence order for non-stiff, medium stiff, and stiff problems. 189 This method is similar to `TSROSWRODASPR`, but satisfies one extra B_{PR} order condition. 190 191 Developer Note: 192 In numerical experiments with ts/tutorials/ex22.c, I (Jed) find this to produce surprisingly poor results. 193 Although the coefficients pass basic smoke tests, I'm not confident it was tabulated correctly in the paper. 194 It would be informative if someone could reproduce tests from the paper and/or reach out to the author to understand why it fails on this test problem. 195 If the method is implemented correctly, doing so might shed light on an additional analysis lens (or further conditions) for robustness on such problems. 196 197 Level: intermediate 198 199 .seealso: [](ch_ts), `TSROSW`, `TSROSWRODASPR` 200 M*/ 201 202 /*MC 203 TSROSWSANDU3 - Three stage third order L-stable Rosenbrock scheme {cite}`sandu_1997` 204 205 By default, the Jacobian is only recomputed once per step. 206 207 The third order method is L-stable, but not stiffly accurate. 208 The second order embedded method is strongly A-stable with R(infty) = 0.5. 209 The internal stages are L-stable. 210 This method is called ROS3 in {cite}`sandu_1997`. 211 212 Level: intermediate 213 214 .seealso: [](ch_ts), `TSROSW`, `TSROSWRODAS3` 215 M*/ 216 217 /*MC 218 TSROSWASSP3P3S1C - A-stable Rosenbrock-W method with SSP explicit part, third order, three stages 219 220 By default, the Jacobian is only recomputed once per step. 221 222 A-stable SPP explicit order 3, 3 stages, CFL 1 (eff = 1/3) 223 224 Level: intermediate 225 226 .seealso: [](ch_ts), `TSROSW`, `TSROSWLASSP3P4S2C`, `TSROSWLLSSP3P4S2C`, `SSP` 227 M*/ 228 229 /*MC 230 TSROSWLASSP3P4S2C - L-stable Rosenbrock-W method with SSP explicit part, third order, four stages 231 232 By default, the Jacobian is only recomputed once per step. 233 234 L-stable (A-stable embedded) SPP explicit order 3, 4 stages, CFL 2 (eff = 1/2) 235 236 Level: intermediate 237 238 .seealso: [](ch_ts), `TSROSW`, `TSROSWASSP3P3S1C`, `TSROSWLLSSP3P4S2C`, `TSSSP` 239 M*/ 240 241 /*MC 242 TSROSWLLSSP3P4S2C - L-stable Rosenbrock-W method with SSP explicit part, third order, four stages 243 244 By default, the Jacobian is only recomputed once per step. 245 246 L-stable (L-stable embedded) SPP explicit order 3, 4 stages, CFL 2 (eff = 1/2) 247 248 Level: intermediate 249 250 .seealso: [](ch_ts), `TSROSW`, `TSROSWASSP3P3S1C`, `TSROSWLASSP3P4S2C`, `TSSSP` 251 M*/ 252 253 /*MC 254 TSROSWGRK4T - four stage, fourth order Rosenbrock (not W) method from Kaps and Rentrop {cite}`kaps1979generalized` 255 256 By default, the Jacobian is only recomputed once per step. 257 258 A(89.3 degrees)-stable, |R(infty)| = 0.454. 259 260 This method does not provide a dense output formula. 261 262 Level: intermediate 263 264 Note: 265 See Section 4 Table 7.2 in {cite}`wanner1996solving` 266 267 .seealso: [](ch_ts), `TSROSW`, `TSROSWSHAMP4`, `TSROSWVELDD4`, `TSROSW4L` 268 M*/ 269 270 /*MC 271 TSROSWSHAMP4 - four stage, fourth order Rosenbrock (not W) method from Shampine {cite}`shampine1982implementation` 272 273 By default, the Jacobian is only recomputed once per step. 274 275 A-stable, |R(infty)| = 1/3. 276 277 This method does not provide a dense output formula. 278 279 Level: intermediate 280 281 Note: 282 See Section 4 Table 7.2 in {cite}`wanner1996solving` 283 284 .seealso: [](ch_ts), `TSROSW`, `TSROSWGRK4T`, `TSROSWVELDD4`, `TSROSW4L` 285 M*/ 286 287 /*MC 288 TSROSWVELDD4 - four stage, fourth order Rosenbrock (not W) method from van Veldhuizen {cite}`veldhuizen1984d` 289 290 By default, the Jacobian is only recomputed once per step. 291 292 A(89.5 degrees)-stable, |R(infty)| = 0.24. 293 294 This method does not provide a dense output formula. 295 296 Level: intermediate 297 298 Note: 299 See Section 4 Table 7.2 in {cite}`wanner1996solving` 300 301 .seealso: [](ch_ts), `TSROSW`, `TSROSWGRK4T`, `TSROSWSHAMP4`, `TSROSW4L` 302 M*/ 303 304 /*MC 305 TSROSW4L - four stage, fourth order Rosenbrock (not W) method 306 307 By default, the Jacobian is only recomputed once per step. 308 309 A-stable and L-stable 310 311 This method does not provide a dense output formula. 312 313 Level: intermediate 314 315 Note: 316 See Section 4 Table 7.2 in {cite}`wanner1996solving` 317 318 .seealso: [](ch_ts), `TSROSW`, `TSROSWGRK4T`, `TSROSWSHAMP4`, `TSROSW4L` 319 M*/ 320 321 /*@C 322 TSRosWRegisterAll - Registers all of the Rosenbrock-W methods in `TSROSW` 323 324 Not Collective, but should be called by all MPI processes which will need the schemes to be registered 325 326 Level: advanced 327 328 .seealso: [](ch_ts), `TSROSW`, `TSRosWRegisterDestroy()` 329 @*/ 330 PetscErrorCode TSRosWRegisterAll(void) 331 { 332 PetscFunctionBegin; 333 if (TSRosWRegisterAllCalled) PetscFunctionReturn(PETSC_SUCCESS); 334 TSRosWRegisterAllCalled = PETSC_TRUE; 335 336 { 337 const PetscReal A = 0; 338 const PetscReal Gamma = 1; 339 const PetscReal b = 1; 340 const PetscReal binterpt = 1; 341 342 PetscCall(TSRosWRegister(TSROSWTHETA1, 1, 1, &A, &Gamma, &b, NULL, 1, &binterpt)); 343 } 344 345 { 346 const PetscReal A = 0; 347 const PetscReal Gamma = 0.5; 348 const PetscReal b = 1; 349 const PetscReal binterpt = 1; 350 351 PetscCall(TSRosWRegister(TSROSWTHETA2, 2, 1, &A, &Gamma, &b, NULL, 1, &binterpt)); 352 } 353 354 { 355 /*const PetscReal g = 1. + 1./PetscSqrtReal(2.0); Direct evaluation: 1.707106781186547524401. Used for setting up arrays of values known at compile time below. */ 356 const PetscReal A[2][2] = { 357 {0, 0}, 358 {1., 0} 359 }; 360 const PetscReal Gamma[2][2] = { 361 {1.707106781186547524401, 0 }, 362 {-2. * 1.707106781186547524401, 1.707106781186547524401} 363 }; 364 const PetscReal b[2] = {0.5, 0.5}; 365 const PetscReal b1[2] = {1.0, 0.0}; 366 PetscReal binterpt[2][2]; 367 binterpt[0][0] = 1.707106781186547524401 - 1.0; 368 binterpt[1][0] = 2.0 - 1.707106781186547524401; 369 binterpt[0][1] = 1.707106781186547524401 - 1.5; 370 binterpt[1][1] = 1.5 - 1.707106781186547524401; 371 372 PetscCall(TSRosWRegister(TSROSW2P, 2, 2, &A[0][0], &Gamma[0][0], b, b1, 2, &binterpt[0][0])); 373 } 374 { 375 /*const PetscReal g = 1. - 1./PetscSqrtReal(2.0); Direct evaluation: 0.2928932188134524755992. Used for setting up arrays of values known at compile time below. */ 376 const PetscReal A[2][2] = { 377 {0, 0}, 378 {1., 0} 379 }; 380 const PetscReal Gamma[2][2] = { 381 {0.2928932188134524755992, 0 }, 382 {-2. * 0.2928932188134524755992, 0.2928932188134524755992} 383 }; 384 const PetscReal b[2] = {0.5, 0.5}; 385 const PetscReal b1[2] = {1.0, 0.0}; 386 PetscReal binterpt[2][2]; 387 binterpt[0][0] = 0.2928932188134524755992 - 1.0; 388 binterpt[1][0] = 2.0 - 0.2928932188134524755992; 389 binterpt[0][1] = 0.2928932188134524755992 - 1.5; 390 binterpt[1][1] = 1.5 - 0.2928932188134524755992; 391 392 PetscCall(TSRosWRegister(TSROSW2M, 2, 2, &A[0][0], &Gamma[0][0], b, b1, 2, &binterpt[0][0])); 393 } 394 { 395 /*const PetscReal g = 7.8867513459481287e-01; Directly written in-place below */ 396 PetscReal binterpt[3][2]; 397 const PetscReal A[3][3] = { 398 {0, 0, 0}, 399 {1.5773502691896257e+00, 0, 0}, 400 {0.5, 0, 0} 401 }; 402 const PetscReal Gamma[3][3] = { 403 {7.8867513459481287e-01, 0, 0 }, 404 {-1.5773502691896257e+00, 7.8867513459481287e-01, 0 }, 405 {-6.7075317547305480e-01, -1.7075317547305482e-01, 7.8867513459481287e-01} 406 }; 407 const PetscReal b[3] = {1.0566243270259355e-01, 4.9038105676657971e-02, 8.4529946162074843e-01}; 408 const PetscReal b2[3] = {-1.7863279495408180e-01, 1. / 3., 8.4529946162074843e-01}; 409 410 binterpt[0][0] = -0.8094010767585034; 411 binterpt[1][0] = -0.5; 412 binterpt[2][0] = 2.3094010767585034; 413 binterpt[0][1] = 0.9641016151377548; 414 binterpt[1][1] = 0.5; 415 binterpt[2][1] = -1.4641016151377548; 416 417 PetscCall(TSRosWRegister(TSROSWRA3PW, 3, 3, &A[0][0], &Gamma[0][0], b, b2, 2, &binterpt[0][0])); 418 } 419 { 420 PetscReal binterpt[4][3]; 421 /*const PetscReal g = 4.3586652150845900e-01; Directly written in-place below */ 422 const PetscReal A[4][4] = { 423 {0, 0, 0, 0}, 424 {8.7173304301691801e-01, 0, 0, 0}, 425 {8.4457060015369423e-01, -1.1299064236484185e-01, 0, 0}, 426 {0, 0, 1., 0} 427 }; 428 const PetscReal Gamma[4][4] = { 429 {4.3586652150845900e-01, 0, 0, 0 }, 430 {-8.7173304301691801e-01, 4.3586652150845900e-01, 0, 0 }, 431 {-9.0338057013044082e-01, 5.4180672388095326e-02, 4.3586652150845900e-01, 0 }, 432 {2.4212380706095346e-01, -1.2232505839045147e+00, 5.4526025533510214e-01, 4.3586652150845900e-01} 433 }; 434 const PetscReal b[4] = {2.4212380706095346e-01, -1.2232505839045147e+00, 1.5452602553351020e+00, 4.3586652150845900e-01}; 435 const PetscReal b2[4] = {3.7810903145819369e-01, -9.6042292212423178e-02, 5.0000000000000000e-01, 2.1793326075422950e-01}; 436 437 binterpt[0][0] = 1.0564298455794094; 438 binterpt[1][0] = 2.296429974281067; 439 binterpt[2][0] = -1.307599564525376; 440 binterpt[3][0] = -1.045260255335102; 441 binterpt[0][1] = -1.3864882699759573; 442 binterpt[1][1] = -8.262611700275677; 443 binterpt[2][1] = 7.250979895056055; 444 binterpt[3][1] = 2.398120075195581; 445 binterpt[0][2] = 0.5721822314575016; 446 binterpt[1][2] = 4.742931142090097; 447 binterpt[2][2] = -4.398120075195578; 448 binterpt[3][2] = -0.9169932983520199; 449 450 PetscCall(TSRosWRegister(TSROSWRA34PW2, 3, 4, &A[0][0], &Gamma[0][0], b, b2, 3, &binterpt[0][0])); 451 } 452 { 453 /* const PetscReal g = 4.3586652150845900e-01; Directly written in-place below */ 454 const PetscReal A[4][4] = { 455 {0, 0, 0, 0}, 456 {8.7173304301691801e-01, 0, 0, 0}, 457 {1.4722022879435914e+00, -3.1840250568090289e-01, 0, 0}, 458 {8.1505192016694938e-01, 5.0000000000000000e-01, -3.1505192016694938e-01, 0} 459 }; 460 const PetscReal Gamma[4][4] = { 461 {4.3586652150845900e-01, 0, 0, 0 }, 462 {-8.7173304301691801e-01, 4.3586652150845900e-01, 0, 0 }, 463 {-1.2855347382089872e+00, 5.0507005541550687e-01, 4.3586652150845900e-01, 0 }, 464 {-4.8201449182864348e-01, 2.1793326075422950e-01, -1.7178529043404503e-01, 4.3586652150845900e-01} 465 }; 466 const PetscReal b[4] = {3.3303742833830591e-01, 7.1793326075422947e-01, -4.8683721060099439e-01, 4.3586652150845900e-01}; 467 const PetscReal b2[4] = {2.5000000000000000e-01, 7.4276119608319180e-01, -3.1472922970066219e-01, 3.2196803361747034e-01}; 468 469 PetscCall(TSRosWRegister(TSROSWR34PRW, 3, 4, &A[0][0], &Gamma[0][0], b, b2, 0, NULL)); 470 } 471 { 472 /* const PetscReal g = 4.3586652150845900e-01; Directly written in-place below */ 473 const PetscReal A[4][4] = { 474 {0, 0, 0, 0}, 475 {1.3075995645253771e+00, 0, 0, 0}, 476 {0.5, 0.5, 0, 0}, 477 {0.5, 0.5, 0, 0} 478 }; 479 const PetscReal Gamma[4][4] = { 480 {4.3586652150845900e-01, 0, 0, 0 }, 481 {-1.3075995645253771e+00, 4.3586652150845900e-01, 0, 0 }, 482 {-7.0988575860972170e-01, -5.5996735960277766e-01, 4.3586652150845900e-01, 0 }, 483 {-1.5550856807552085e-01, -9.5388516575112225e-01, 6.7352721231818413e-01, 4.3586652150845900e-01} 484 }; 485 const PetscReal b[4] = {3.4449143192447917e-01, -4.5388516575112231e-01, 6.7352721231818413e-01, 4.3586652150845900e-01}; 486 const PetscReal b2[4] = {5.0000000000000000e-01, -2.5738812086522078e-01, 4.3542008724775044e-01, 3.2196803361747034e-01}; 487 488 PetscCall(TSRosWRegister(TSROSWR3PRL2, 3, 4, &A[0][0], &Gamma[0][0], b, b2, 0, NULL)); 489 } 490 { 491 /* const PetscReal g = 0.5; Directly written in-place below */ 492 const PetscReal A[4][4] = { 493 {0, 0, 0, 0}, 494 {0, 0, 0, 0}, 495 {1., 0, 0, 0}, 496 {0.75, -0.25, 0.5, 0} 497 }; 498 const PetscReal Gamma[4][4] = { 499 {0.5, 0, 0, 0 }, 500 {1., 0.5, 0, 0 }, 501 {-0.25, -0.25, 0.5, 0 }, 502 {1. / 12, 1. / 12, -2. / 3, 0.5} 503 }; 504 const PetscReal b[4] = {5. / 6, -1. / 6, -1. / 6, 0.5}; 505 const PetscReal b2[4] = {0.75, -0.25, 0.5, 0}; 506 507 PetscCall(TSRosWRegister(TSROSWRODAS3, 3, 4, &A[0][0], &Gamma[0][0], b, b2, 0, NULL)); 508 } 509 { 510 /* const PetscReal g = 0.25; Directly written in-place below */ 511 const PetscReal A[6][6] = { 512 {0, 0, 0, 0, 0, 0}, 513 {0.75, 0, 0, 0, 0, 0}, 514 {7.5162877593868457e-02, 2.4837122406131545e-02, 0, 0, 0, 0}, 515 {1.6532708886396510e+00, 2.1545706385445562e-01, -1.3157488872766792e+00, 0, 0, 0}, 516 {1.9385003738039885e+01, 1.2007117225835324e+00, -1.9337924059522791e+01, -2.4779140110062559e-01, 0, 0}, 517 {-7.3844531665375115e+00, -3.0593419030174646e-01, 7.8622074209377981e+00, 5.7817993590145966e-01, 0.25, 0} 518 }; 519 const PetscReal Gamma[6][6] = { 520 {0.25, 0, 0, 0, 0, 0 }, 521 {-7.5000000000000000e-01, 0.25, 0, 0, 0, 0 }, 522 {-8.8644359075349941e-02, -2.8688974257983398e-02, 0.25, 0, 0, 0 }, 523 {-4.8470034585330284e+00, -3.1583244269672095e-01, 4.9536568360123221e+00, 0.25, 0, 0 }, 524 {-2.6769456904577400e+01, -1.5066459128852787e+00, 2.7200131480460591e+01, 8.2597133700208525e-01, 0.25, 0 }, 525 {6.5876206496361416e+00, 3.6807059172993878e-01, -6.7423520694658121e+00, -1.0619631475741095e-01, -3.5714285714285715e-01, 0.25} 526 }; 527 const PetscReal b[6] = {-7.9683251690137014e-01, 6.2136401428192344e-02, 1.1198553514719862e+00, 4.7198362114404874e-01, -1.0714285714285714e-01, 0.25}; 528 const PetscReal b2[6] = {-7.3844531665375115e+00, -3.0593419030174646e-01, 7.8622074209377981e+00, 5.7817993590145966e-01, 0.25, 0.0}; 529 530 PetscCall(TSRosWRegister(TSROSWRODASPR, 4, 6, &A[0][0], &Gamma[0][0], b, b2, 0, NULL)); 531 } 532 { 533 /* const PetscReal g = 0.3125; Directly written in-place below */ 534 const PetscReal A[6][6] = { 535 {0, 0, 0, 0, 0, 0}, 536 {9.3750000000000000e-01, 0, 0, 0, 0, 0}, 537 {-4.7145892646261345e-02, 5.4531286650471122e-01, 0, 0, 0, 0}, 538 {4.6915543899742240e-01, 4.4490537602383673e-01, -2.2498239334061121e-01, 0, 0, 0}, 539 {1.0950372887345903e+00, 6.3223023457294381e-01, -8.9232966090485821e-01, 1.6506213759732410e-01, 0, 0}, 540 {-1.7746585073632790e-01, -5.8241418952602364e-01, 6.8180612588238165e-01, 7.6557391437996980e-01, 3.1250000000000000e-01, 0} 541 }; 542 const PetscReal Gamma[6][6] = { 543 {0.3125, 0, 0, 0, 0, 0 }, 544 {-9.3750000000000000e-01, 0.3125, 0, 0, 0, 0 }, 545 {-9.7580572085994507e-02, -5.8666328499964138e-01, 0.3125, 0, 0, 0 }, 546 {-4.9407065013256957e-01, -5.6819726428975503e-01, 5.0318949274167679e-01, 0.3125, 0, 0 }, 547 {-1.2725031394709183e+00, -1.2146444240989676e+00, 1.5741357867872399e+00, 6.0051177678264578e-01, 0.3125, 0 }, 548 {6.9690744901421153e-01, 6.2237005730756434e-01, -1.1553701989197045e+00, 1.8350029013386296e-01, -6.5990759753593431e-01, 0.3125} 549 }; 550 const PetscReal b[6] = {5.1944159827788361e-01, 3.9955867781540699e-02, -4.7356407303732290e-01, 9.4907420451383284e-01, -3.4740759753593431e-01, 0.3125}; 551 const PetscReal b2[6] = {-1.7746585073632790e-01, -5.8241418952602364e-01, 6.8180612588238165e-01, 7.6557391437996980e-01, 0.3125, 0}; 552 553 PetscCall(TSRosWRegister(TSROSWRODASPR2, 4, 6, &A[0][0], &Gamma[0][0], b, b2, 0, NULL)); 554 } 555 { 556 /*const PetscReal g = 0.43586652150845899941601945119356; Directly written in-place below */ 557 const PetscReal A[3][3] = { 558 {0, 0, 0}, 559 {0.43586652150845899941601945119356, 0, 0}, 560 {0.43586652150845899941601945119356, 0, 0} 561 }; 562 const PetscReal Gamma[3][3] = { 563 {0.43586652150845899941601945119356, 0, 0 }, 564 {-0.19294655696029095575009695436041, 0.43586652150845899941601945119356, 0 }, 565 {0, 1.74927148125794685173529749738960, 0.43586652150845899941601945119356} 566 }; 567 const PetscReal b[3] = {-0.75457412385404315829818998646589, 1.94100407061964420292840123379419, -0.18642994676560104463021124732829}; 568 const PetscReal b2[3] = {-1.53358745784149585370766523913002, 2.81745131148625772213931745457622, -0.28386385364476186843165221544619}; 569 570 PetscReal binterpt[3][2]; 571 binterpt[0][0] = 3.793692883777660870425141387941; 572 binterpt[1][0] = -2.918692883777660870425141387941; 573 binterpt[2][0] = 0.125; 574 binterpt[0][1] = -0.725741064379812106687651020584; 575 binterpt[1][1] = 0.559074397713145440020984353917; 576 binterpt[2][1] = 0.16666666666666666666666666666667; 577 578 PetscCall(TSRosWRegister(TSROSWSANDU3, 3, 3, &A[0][0], &Gamma[0][0], b, b2, 2, &binterpt[0][0])); 579 } 580 { 581 /*const PetscReal s3 = PetscSqrtReal(3.),g = (3.0+s3)/6.0; 582 * Direct evaluation: s3 = 1.732050807568877293527; 583 * g = 0.7886751345948128822546; 584 * Values are directly inserted below to ensure availability at compile time (compiler warnings otherwise...) */ 585 const PetscReal A[3][3] = { 586 {0, 0, 0}, 587 {1, 0, 0}, 588 {0.25, 0.25, 0} 589 }; 590 const PetscReal Gamma[3][3] = { 591 {0, 0, 0 }, 592 {(-3.0 - 1.732050807568877293527) / 6.0, 0.7886751345948128822546, 0 }, 593 {(-3.0 - 1.732050807568877293527) / 24.0, (-3.0 - 1.732050807568877293527) / 8.0, 0.7886751345948128822546} 594 }; 595 const PetscReal b[3] = {1. / 6., 1. / 6., 2. / 3.}; 596 const PetscReal b2[3] = {1. / 4., 1. / 4., 1. / 2.}; 597 PetscReal binterpt[3][2]; 598 599 binterpt[0][0] = 0.089316397477040902157517886164709; 600 binterpt[1][0] = -0.91068360252295909784248211383529; 601 binterpt[2][0] = 1.8213672050459181956849642276706; 602 binterpt[0][1] = 0.077350269189625764509148780501957; 603 binterpt[1][1] = 1.077350269189625764509148780502; 604 binterpt[2][1] = -1.1547005383792515290182975610039; 605 606 PetscCall(TSRosWRegister(TSROSWASSP3P3S1C, 3, 3, &A[0][0], &Gamma[0][0], b, b2, 2, &binterpt[0][0])); 607 } 608 609 { 610 const PetscReal A[4][4] = { 611 {0, 0, 0, 0}, 612 {1. / 2., 0, 0, 0}, 613 {1. / 2., 1. / 2., 0, 0}, 614 {1. / 6., 1. / 6., 1. / 6., 0} 615 }; 616 const PetscReal Gamma[4][4] = { 617 {1. / 2., 0, 0, 0}, 618 {0.0, 1. / 4., 0, 0}, 619 {-2., -2. / 3., 2. / 3., 0}, 620 {1. / 2., 5. / 36., -2. / 9, 0} 621 }; 622 const PetscReal b[4] = {1. / 6., 1. / 6., 1. / 6., 1. / 2.}; 623 const PetscReal b2[4] = {1. / 8., 3. / 4., 1. / 8., 0}; 624 PetscReal binterpt[4][3]; 625 626 binterpt[0][0] = 6.25; 627 binterpt[1][0] = -30.25; 628 binterpt[2][0] = 1.75; 629 binterpt[3][0] = 23.25; 630 binterpt[0][1] = -9.75; 631 binterpt[1][1] = 58.75; 632 binterpt[2][1] = -3.25; 633 binterpt[3][1] = -45.75; 634 binterpt[0][2] = 3.6666666666666666666666666666667; 635 binterpt[1][2] = -28.333333333333333333333333333333; 636 binterpt[2][2] = 1.6666666666666666666666666666667; 637 binterpt[3][2] = 23.; 638 639 PetscCall(TSRosWRegister(TSROSWLASSP3P4S2C, 3, 4, &A[0][0], &Gamma[0][0], b, b2, 3, &binterpt[0][0])); 640 } 641 642 { 643 const PetscReal A[4][4] = { 644 {0, 0, 0, 0}, 645 {1. / 2., 0, 0, 0}, 646 {1. / 2., 1. / 2., 0, 0}, 647 {1. / 6., 1. / 6., 1. / 6., 0} 648 }; 649 const PetscReal Gamma[4][4] = { 650 {1. / 2., 0, 0, 0}, 651 {0.0, 3. / 4., 0, 0}, 652 {-2. / 3., -23. / 9., 2. / 9., 0}, 653 {1. / 18., 65. / 108., -2. / 27, 0} 654 }; 655 const PetscReal b[4] = {1. / 6., 1. / 6., 1. / 6., 1. / 2.}; 656 const PetscReal b2[4] = {3. / 16., 10. / 16., 3. / 16., 0}; 657 PetscReal binterpt[4][3]; 658 659 binterpt[0][0] = 1.6911764705882352941176470588235; 660 binterpt[1][0] = 3.6813725490196078431372549019608; 661 binterpt[2][0] = 0.23039215686274509803921568627451; 662 binterpt[3][0] = -4.6029411764705882352941176470588; 663 binterpt[0][1] = -0.95588235294117647058823529411765; 664 binterpt[1][1] = -6.2401960784313725490196078431373; 665 binterpt[2][1] = -0.31862745098039215686274509803922; 666 binterpt[3][1] = 7.5147058823529411764705882352941; 667 binterpt[0][2] = -0.56862745098039215686274509803922; 668 binterpt[1][2] = 2.7254901960784313725490196078431; 669 binterpt[2][2] = 0.25490196078431372549019607843137; 670 binterpt[3][2] = -2.4117647058823529411764705882353; 671 672 PetscCall(TSRosWRegister(TSROSWLLSSP3P4S2C, 3, 4, &A[0][0], &Gamma[0][0], b, b2, 3, &binterpt[0][0])); 673 } 674 675 { 676 PetscReal A[4][4], Gamma[4][4], b[4], b2[4]; 677 PetscReal binterpt[4][3]; 678 679 Gamma[0][0] = 0.4358665215084589994160194475295062513822671686978816; 680 Gamma[0][1] = 0; 681 Gamma[0][2] = 0; 682 Gamma[0][3] = 0; 683 Gamma[1][0] = -1.997527830934941248426324674704153457289527280554476; 684 Gamma[1][1] = 0.4358665215084589994160194475295062513822671686978816; 685 Gamma[1][2] = 0; 686 Gamma[1][3] = 0; 687 Gamma[2][0] = -1.007948511795029620852002345345404191008352770119903; 688 Gamma[2][1] = -0.004648958462629345562774289390054679806993396798458131; 689 Gamma[2][2] = 0.4358665215084589994160194475295062513822671686978816; 690 Gamma[2][3] = 0; 691 Gamma[3][0] = -0.6685429734233467180451604600279552604364311322650783; 692 Gamma[3][1] = 0.6056625986449338476089525334450053439525178740492984; 693 Gamma[3][2] = -0.9717899277217721234705114616271378792182450260943198; 694 Gamma[3][3] = 0; 695 696 A[0][0] = 0; 697 A[0][1] = 0; 698 A[0][2] = 0; 699 A[0][3] = 0; 700 A[1][0] = 0.8717330430169179988320388950590125027645343373957631; 701 A[1][1] = 0; 702 A[1][2] = 0; 703 A[1][3] = 0; 704 A[2][0] = 0.5275890119763004115618079766722914408876108660811028; 705 A[2][1] = 0.07241098802369958843819203208518599088698057726988732; 706 A[2][2] = 0; 707 A[2][3] = 0; 708 A[3][0] = 0.3990960076760701320627260685975778145384666450351314; 709 A[3][1] = -0.4375576546135194437228463747348862825846903771419953; 710 A[3][2] = 1.038461646937449311660120300601880176655352737312713; 711 A[3][3] = 0; 712 713 b[0] = 0.1876410243467238251612921333138006734899663569186926; 714 b[1] = -0.5952974735769549480478230473706443582188442040780541; 715 b[2] = 0.9717899277217721234705114616271378792182450260943198; 716 b[3] = 0.4358665215084589994160194475295062513822671686978816; 717 718 b2[0] = 0.2147402862233891404862383521089097657790734483804460; 719 b2[1] = -0.4851622638849390928209050538171743017757490232519684; 720 b2[2] = 0.8687250025203875511662123688667549217531982787600080; 721 b2[3] = 0.4016969751411624011684543450940068201770721128357014; 722 723 binterpt[0][0] = 2.2565812720167954547104627844105; 724 binterpt[1][0] = 1.349166413351089573796243820819; 725 binterpt[2][0] = -2.4695174540533503758652847586647; 726 binterpt[3][0] = -0.13623023131453465264142184656474; 727 binterpt[0][1] = -3.0826699111559187902922463354557; 728 binterpt[1][1] = -2.4689115685996042534544925650515; 729 binterpt[2][1] = 5.7428279814696677152129332773553; 730 binterpt[3][1] = -0.19124650171414467146619437684812; 731 binterpt[0][2] = 1.0137296634858471607430756831148; 732 binterpt[1][2] = 0.52444768167155973161042570784064; 733 binterpt[2][2] = -2.3015205996945452158771370439586; 734 binterpt[3][2] = 0.76334325453713832352363565300308; 735 736 PetscCall(TSRosWRegister(TSROSWARK3, 3, 4, &A[0][0], &Gamma[0][0], b, b2, 3, &binterpt[0][0])); 737 } 738 PetscCall(TSRosWRegisterRos4(TSROSWGRK4T, 0.231, PETSC_CURRENT, PETSC_CURRENT, 0, -0.1282612945269037e+01)); 739 PetscCall(TSRosWRegisterRos4(TSROSWSHAMP4, 0.5, PETSC_CURRENT, PETSC_CURRENT, 0, 125. / 108.)); 740 PetscCall(TSRosWRegisterRos4(TSROSWVELDD4, 0.22570811482256823492, PETSC_CURRENT, PETSC_CURRENT, 0, -1.355958941201148)); 741 PetscCall(TSRosWRegisterRos4(TSROSW4L, 0.57282, PETSC_CURRENT, PETSC_CURRENT, 0, -1.093502252409163)); 742 PetscFunctionReturn(PETSC_SUCCESS); 743 } 744 745 /*@C 746 TSRosWRegisterDestroy - Frees the list of schemes that were registered by `TSRosWRegister()`. 747 748 Not Collective 749 750 Level: advanced 751 752 .seealso: [](ch_ts), `TSRosWRegister()`, `TSRosWRegisterAll()` 753 @*/ 754 PetscErrorCode TSRosWRegisterDestroy(void) 755 { 756 RosWTableauLink link; 757 758 PetscFunctionBegin; 759 while ((link = RosWTableauList)) { 760 RosWTableau t = &link->tab; 761 RosWTableauList = link->next; 762 PetscCall(PetscFree5(t->A, t->Gamma, t->b, t->ASum, t->GammaSum)); 763 PetscCall(PetscFree5(t->At, t->bt, t->GammaInv, t->GammaZeroDiag, t->GammaExplicitCorr)); 764 PetscCall(PetscFree2(t->bembed, t->bembedt)); 765 PetscCall(PetscFree(t->binterpt)); 766 PetscCall(PetscFree(t->name)); 767 PetscCall(PetscFree(link)); 768 } 769 TSRosWRegisterAllCalled = PETSC_FALSE; 770 PetscFunctionReturn(PETSC_SUCCESS); 771 } 772 773 /*@C 774 TSRosWInitializePackage - This function initializes everything in the `TSROSW` package. It is called 775 from `TSInitializePackage()`. 776 777 Level: developer 778 779 .seealso: [](ch_ts), `TSROSW`, `PetscInitialize()`, `TSRosWFinalizePackage()` 780 @*/ 781 PetscErrorCode TSRosWInitializePackage(void) 782 { 783 PetscFunctionBegin; 784 if (TSRosWPackageInitialized) PetscFunctionReturn(PETSC_SUCCESS); 785 TSRosWPackageInitialized = PETSC_TRUE; 786 PetscCall(TSRosWRegisterAll()); 787 PetscCall(PetscRegisterFinalize(TSRosWFinalizePackage)); 788 PetscFunctionReturn(PETSC_SUCCESS); 789 } 790 791 /*@C 792 TSRosWFinalizePackage - This function destroys everything in the `TSROSW` package. It is 793 called from `PetscFinalize()`. 794 795 Level: developer 796 797 .seealso: [](ch_ts), `TSROSW`, `PetscFinalize()`, `TSRosWInitializePackage()` 798 @*/ 799 PetscErrorCode TSRosWFinalizePackage(void) 800 { 801 PetscFunctionBegin; 802 TSRosWPackageInitialized = PETSC_FALSE; 803 PetscCall(TSRosWRegisterDestroy()); 804 PetscFunctionReturn(PETSC_SUCCESS); 805 } 806 807 /*@C 808 TSRosWRegister - register a `TSROSW`, Rosenbrock W scheme by providing the entries in the Butcher tableau and optionally embedded approximations and interpolation 809 810 Not Collective, but the same schemes should be registered on all processes on which they will be used 811 812 Input Parameters: 813 + name - identifier for method 814 . order - approximation order of method 815 . s - number of stages, this is the dimension of the matrices below 816 . A - Table of propagated stage coefficients (dimension s*s, row-major), strictly lower triangular 817 . Gamma - Table of coefficients in implicit stage equations (dimension s*s, row-major), lower triangular with nonzero diagonal 818 . b - Step completion table (dimension s) 819 . bembed - Step completion table for a scheme of order one less (dimension s, NULL if no embedded scheme is available) 820 . pinterp - Order of the interpolation scheme, equal to the number of columns of binterpt 821 - binterpt - Coefficients of the interpolation formula (dimension s*pinterp) 822 823 Level: advanced 824 825 Note: 826 Several Rosenbrock W methods are provided, this function is only needed to create new methods. 827 828 .seealso: [](ch_ts), `TSROSW` 829 @*/ 830 PetscErrorCode TSRosWRegister(TSRosWType name, PetscInt order, PetscInt s, const PetscReal A[], const PetscReal Gamma[], const PetscReal b[], const PetscReal bembed[], PetscInt pinterp, const PetscReal binterpt[]) 831 { 832 RosWTableauLink link; 833 RosWTableau t; 834 PetscInt i, j, k; 835 PetscScalar *GammaInv; 836 837 PetscFunctionBegin; 838 PetscAssertPointer(name, 1); 839 PetscAssertPointer(A, 4); 840 PetscAssertPointer(Gamma, 5); 841 PetscAssertPointer(b, 6); 842 if (bembed) PetscAssertPointer(bembed, 7); 843 844 PetscCall(TSRosWInitializePackage()); 845 PetscCall(PetscNew(&link)); 846 t = &link->tab; 847 PetscCall(PetscStrallocpy(name, &t->name)); 848 t->order = order; 849 t->s = s; 850 PetscCall(PetscMalloc5(s * s, &t->A, s * s, &t->Gamma, s, &t->b, s, &t->ASum, s, &t->GammaSum)); 851 PetscCall(PetscMalloc5(s * s, &t->At, s, &t->bt, s * s, &t->GammaInv, s, &t->GammaZeroDiag, s * s, &t->GammaExplicitCorr)); 852 PetscCall(PetscArraycpy(t->A, A, s * s)); 853 PetscCall(PetscArraycpy(t->Gamma, Gamma, s * s)); 854 PetscCall(PetscArraycpy(t->GammaExplicitCorr, Gamma, s * s)); 855 PetscCall(PetscArraycpy(t->b, b, s)); 856 if (bembed) { 857 PetscCall(PetscMalloc2(s, &t->bembed, s, &t->bembedt)); 858 PetscCall(PetscArraycpy(t->bembed, bembed, s)); 859 } 860 for (i = 0; i < s; i++) { 861 t->ASum[i] = 0; 862 t->GammaSum[i] = 0; 863 for (j = 0; j < s; j++) { 864 t->ASum[i] += A[i * s + j]; 865 t->GammaSum[i] += Gamma[i * s + j]; 866 } 867 } 868 PetscCall(PetscMalloc1(s * s, &GammaInv)); /* Need to use Scalar for inverse, then convert back to Real */ 869 for (i = 0; i < s * s; i++) GammaInv[i] = Gamma[i]; 870 for (i = 0; i < s; i++) { 871 if (Gamma[i * s + i] == 0.0) { 872 GammaInv[i * s + i] = 1.0; 873 t->GammaZeroDiag[i] = PETSC_TRUE; 874 } else { 875 t->GammaZeroDiag[i] = PETSC_FALSE; 876 } 877 } 878 879 switch (s) { 880 case 1: 881 GammaInv[0] = 1. / GammaInv[0]; 882 break; 883 case 2: 884 PetscCall(PetscKernel_A_gets_inverse_A_2(GammaInv, 0, PETSC_FALSE, NULL)); 885 break; 886 case 3: 887 PetscCall(PetscKernel_A_gets_inverse_A_3(GammaInv, 0, PETSC_FALSE, NULL)); 888 break; 889 case 4: 890 PetscCall(PetscKernel_A_gets_inverse_A_4(GammaInv, 0, PETSC_FALSE, NULL)); 891 break; 892 case 5: { 893 PetscInt ipvt5[5]; 894 MatScalar work5[5 * 5]; 895 PetscCall(PetscKernel_A_gets_inverse_A_5(GammaInv, ipvt5, work5, 0, PETSC_FALSE, NULL)); 896 break; 897 } 898 case 6: 899 PetscCall(PetscKernel_A_gets_inverse_A_6(GammaInv, 0, PETSC_FALSE, NULL)); 900 break; 901 case 7: 902 PetscCall(PetscKernel_A_gets_inverse_A_7(GammaInv, 0, PETSC_FALSE, NULL)); 903 break; 904 default: 905 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Not implemented for %" PetscInt_FMT " stages", s); 906 } 907 for (i = 0; i < s * s; i++) t->GammaInv[i] = PetscRealPart(GammaInv[i]); 908 PetscCall(PetscFree(GammaInv)); 909 910 for (i = 0; i < s; i++) { 911 for (k = 0; k < i + 1; k++) { 912 t->GammaExplicitCorr[i * s + k] = (t->GammaExplicitCorr[i * s + k]) * (t->GammaInv[k * s + k]); 913 for (j = k + 1; j < i + 1; j++) t->GammaExplicitCorr[i * s + k] += (t->GammaExplicitCorr[i * s + j]) * (t->GammaInv[j * s + k]); 914 } 915 } 916 917 for (i = 0; i < s; i++) { 918 for (j = 0; j < s; j++) { 919 t->At[i * s + j] = 0; 920 for (k = 0; k < s; k++) t->At[i * s + j] += t->A[i * s + k] * t->GammaInv[k * s + j]; 921 } 922 t->bt[i] = 0; 923 for (j = 0; j < s; j++) t->bt[i] += t->b[j] * t->GammaInv[j * s + i]; 924 if (bembed) { 925 t->bembedt[i] = 0; 926 for (j = 0; j < s; j++) t->bembedt[i] += t->bembed[j] * t->GammaInv[j * s + i]; 927 } 928 } 929 t->ccfl = 1.0; /* Fix this */ 930 931 t->pinterp = pinterp; 932 PetscCall(PetscMalloc1(s * pinterp, &t->binterpt)); 933 PetscCall(PetscArraycpy(t->binterpt, binterpt, s * pinterp)); 934 link->next = RosWTableauList; 935 RosWTableauList = link; 936 PetscFunctionReturn(PETSC_SUCCESS); 937 } 938 939 /*@C 940 TSRosWRegisterRos4 - register a fourth order Rosenbrock scheme by providing parameter choices 941 942 Not Collective, but the same schemes should be registered on all processes on which they will be used 943 944 Input Parameters: 945 + name - identifier for method 946 . gamma - leading coefficient (diagonal entry) 947 . a2 - design parameter, see Table 7.2 of {cite}`wanner1996solving` 948 . a3 - design parameter or `PETSC_DETERMINE` to satisfy one of the order five conditions (Eq 7.22) 949 . b3 - design parameter, see Table 7.2 of {cite}`wanner1996solving` 950 - e4 - design parameter for embedded method, see coefficient E4 in ros4.f code from Hairer 951 952 Level: developer 953 954 Notes: 955 This routine encodes the design of fourth order Rosenbrock methods as described in {cite}`wanner1996solving` 956 It is used here to implement several methods from the book and can be used to experiment with new methods. 957 It was written this way instead of by copying coefficients in order to provide better than double precision satisfaction of the order conditions. 958 959 .seealso: [](ch_ts), `TSROSW`, `TSRosWRegister()` 960 @*/ 961 PetscErrorCode TSRosWRegisterRos4(TSRosWType name, PetscReal gamma, PetscReal a2, PetscReal a3, PetscReal b3, PetscReal e4) 962 { 963 /* Declare numeric constants so they can be quad precision without being truncated at double */ 964 const PetscReal one = 1, two = 2, three = 3, four = 4, five = 5, six = 6, eight = 8, twelve = 12, twenty = 20, twentyfour = 24, p32 = one / six - gamma + gamma * gamma, p42 = one / eight - gamma / three, p43 = one / twelve - gamma / three, p44 = one / twentyfour - gamma / two + three / two * gamma * gamma - gamma * gamma * gamma, p56 = one / twenty - gamma / four; 965 PetscReal a4, a32, a42, a43, b1, b2, b4, beta2p, beta3p, beta4p, beta32, beta42, beta43, beta32beta2p, beta4jbetajp; 966 PetscReal A[4][4], Gamma[4][4], b[4], bm[4]; 967 PetscScalar M[3][3], rhs[3]; 968 969 PetscFunctionBegin; 970 /* Step 1: choose Gamma (input) */ 971 /* Step 2: choose a2,a3,a4; b1,b2,b3,b4 to satisfy order conditions */ 972 if (a3 == (PetscReal)PETSC_DEFAULT || a3 == (PetscReal)PETSC_DETERMINE) a3 = (one / five - a2 / four) / (one / four - a2 / three); /* Eq 7.22 */ 973 a4 = a3; /* consequence of 7.20 */ 974 975 /* Solve order conditions 7.15a, 7.15c, 7.15e */ 976 M[0][0] = one; 977 M[0][1] = one; 978 M[0][2] = one; /* 7.15a */ 979 M[1][0] = 0.0; 980 M[1][1] = a2 * a2; 981 M[1][2] = a4 * a4; /* 7.15c */ 982 M[2][0] = 0.0; 983 M[2][1] = a2 * a2 * a2; 984 M[2][2] = a4 * a4 * a4; /* 7.15e */ 985 rhs[0] = one - b3; 986 rhs[1] = one / three - a3 * a3 * b3; 987 rhs[2] = one / four - a3 * a3 * a3 * b3; 988 PetscCall(PetscKernel_A_gets_inverse_A_3(&M[0][0], 0, PETSC_FALSE, NULL)); 989 b1 = PetscRealPart(M[0][0] * rhs[0] + M[0][1] * rhs[1] + M[0][2] * rhs[2]); 990 b2 = PetscRealPart(M[1][0] * rhs[0] + M[1][1] * rhs[1] + M[1][2] * rhs[2]); 991 b4 = PetscRealPart(M[2][0] * rhs[0] + M[2][1] * rhs[1] + M[2][2] * rhs[2]); 992 993 /* Step 3 */ 994 beta43 = (p56 - a2 * p43) / (b4 * a3 * a3 * (a3 - a2)); /* 7.21 */ 995 beta32beta2p = p44 / (b4 * beta43); /* 7.15h */ 996 beta4jbetajp = (p32 - b3 * beta32beta2p) / b4; 997 M[0][0] = b2; 998 M[0][1] = b3; 999 M[0][2] = b4; 1000 M[1][0] = a4 * a4 * beta32beta2p - a3 * a3 * beta4jbetajp; 1001 M[1][1] = a2 * a2 * beta4jbetajp; 1002 M[1][2] = -a2 * a2 * beta32beta2p; 1003 M[2][0] = b4 * beta43 * a3 * a3 - p43; 1004 M[2][1] = -b4 * beta43 * a2 * a2; 1005 M[2][2] = 0; 1006 rhs[0] = one / two - gamma; 1007 rhs[1] = 0; 1008 rhs[2] = -a2 * a2 * p32; 1009 PetscCall(PetscKernel_A_gets_inverse_A_3(&M[0][0], 0, PETSC_FALSE, NULL)); 1010 beta2p = PetscRealPart(M[0][0] * rhs[0] + M[0][1] * rhs[1] + M[0][2] * rhs[2]); 1011 beta3p = PetscRealPart(M[1][0] * rhs[0] + M[1][1] * rhs[1] + M[1][2] * rhs[2]); 1012 beta4p = PetscRealPart(M[2][0] * rhs[0] + M[2][1] * rhs[1] + M[2][2] * rhs[2]); 1013 1014 /* Step 4: back-substitute */ 1015 beta32 = beta32beta2p / beta2p; 1016 beta42 = (beta4jbetajp - beta43 * beta3p) / beta2p; 1017 1018 /* Step 5: 7.15f and 7.20, then 7.16 */ 1019 a43 = 0; 1020 a32 = p42 / (b3 * a3 * beta2p + b4 * a4 * beta2p); 1021 a42 = a32; 1022 1023 A[0][0] = 0; 1024 A[0][1] = 0; 1025 A[0][2] = 0; 1026 A[0][3] = 0; 1027 A[1][0] = a2; 1028 A[1][1] = 0; 1029 A[1][2] = 0; 1030 A[1][3] = 0; 1031 A[2][0] = a3 - a32; 1032 A[2][1] = a32; 1033 A[2][2] = 0; 1034 A[2][3] = 0; 1035 A[3][0] = a4 - a43 - a42; 1036 A[3][1] = a42; 1037 A[3][2] = a43; 1038 A[3][3] = 0; 1039 Gamma[0][0] = gamma; 1040 Gamma[0][1] = 0; 1041 Gamma[0][2] = 0; 1042 Gamma[0][3] = 0; 1043 Gamma[1][0] = beta2p - A[1][0]; 1044 Gamma[1][1] = gamma; 1045 Gamma[1][2] = 0; 1046 Gamma[1][3] = 0; 1047 Gamma[2][0] = beta3p - beta32 - A[2][0]; 1048 Gamma[2][1] = beta32 - A[2][1]; 1049 Gamma[2][2] = gamma; 1050 Gamma[2][3] = 0; 1051 Gamma[3][0] = beta4p - beta42 - beta43 - A[3][0]; 1052 Gamma[3][1] = beta42 - A[3][1]; 1053 Gamma[3][2] = beta43 - A[3][2]; 1054 Gamma[3][3] = gamma; 1055 b[0] = b1; 1056 b[1] = b2; 1057 b[2] = b3; 1058 b[3] = b4; 1059 1060 /* Construct embedded formula using given e4. We are solving Equation 7.18. */ 1061 bm[3] = b[3] - e4 * gamma; /* using definition of E4 */ 1062 bm[2] = (p32 - beta4jbetajp * bm[3]) / (beta32 * beta2p); /* fourth row of 7.18 */ 1063 bm[1] = (one / two - gamma - beta3p * bm[2] - beta4p * bm[3]) / beta2p; /* second row */ 1064 bm[0] = one - bm[1] - bm[2] - bm[3]; /* first row */ 1065 1066 { 1067 const PetscReal misfit = a2 * a2 * bm[1] + a3 * a3 * bm[2] + a4 * a4 * bm[3] - one / three; 1068 PetscCheck(PetscAbs(misfit) <= PETSC_SMALL, PETSC_COMM_SELF, PETSC_ERR_SUP, "Assumptions violated, could not construct a third order embedded method"); 1069 } 1070 PetscCall(TSRosWRegister(name, 4, 4, &A[0][0], &Gamma[0][0], b, bm, 0, NULL)); 1071 PetscFunctionReturn(PETSC_SUCCESS); 1072 } 1073 1074 /* 1075 The step completion formula is 1076 1077 x1 = x0 + b^T Y 1078 1079 where Y is the multi-vector of stages corrections. This function can be called before or after ts->vec_sol has been 1080 updated. Suppose we have a completion formula b and an embedded formula be of different order. We can write 1081 1082 x1e = x0 + be^T Y 1083 = x1 - b^T Y + be^T Y 1084 = x1 + (be - b)^T Y 1085 1086 so we can evaluate the method of different order even after the step has been optimistically completed. 1087 */ 1088 static PetscErrorCode TSEvaluateStep_RosW(TS ts, PetscInt order, Vec U, PetscBool *done) 1089 { 1090 TS_RosW *ros = (TS_RosW *)ts->data; 1091 RosWTableau tab = ros->tableau; 1092 PetscScalar *w = ros->work; 1093 PetscInt i; 1094 1095 PetscFunctionBegin; 1096 if (order == tab->order) { 1097 if (ros->status == TS_STEP_INCOMPLETE) { /* Use standard completion formula */ 1098 PetscCall(VecCopy(ts->vec_sol, U)); 1099 for (i = 0; i < tab->s; i++) w[i] = tab->bt[i]; 1100 PetscCall(VecMAXPY(U, tab->s, w, ros->Y)); 1101 } else PetscCall(VecCopy(ts->vec_sol, U)); 1102 if (done) *done = PETSC_TRUE; 1103 PetscFunctionReturn(PETSC_SUCCESS); 1104 } else if (order == tab->order - 1) { 1105 if (!tab->bembedt) goto unavailable; 1106 if (ros->status == TS_STEP_INCOMPLETE) { /* Use embedded completion formula */ 1107 PetscCall(VecCopy(ts->vec_sol, U)); 1108 for (i = 0; i < tab->s; i++) w[i] = tab->bembedt[i]; 1109 PetscCall(VecMAXPY(U, tab->s, w, ros->Y)); 1110 } else { /* Use rollback-and-recomplete formula (bembedt - bt) */ 1111 for (i = 0; i < tab->s; i++) w[i] = tab->bembedt[i] - tab->bt[i]; 1112 PetscCall(VecCopy(ts->vec_sol, U)); 1113 PetscCall(VecMAXPY(U, tab->s, w, ros->Y)); 1114 } 1115 if (done) *done = PETSC_TRUE; 1116 PetscFunctionReturn(PETSC_SUCCESS); 1117 } 1118 unavailable: 1119 if (done) *done = PETSC_FALSE; 1120 else 1121 SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "Rosenbrock-W '%s' of order %" PetscInt_FMT " cannot evaluate step at order %" PetscInt_FMT ". Consider using -ts_adapt_type none or a different method that has an embedded estimate.", tab->name, 1122 tab->order, order); 1123 PetscFunctionReturn(PETSC_SUCCESS); 1124 } 1125 1126 static PetscErrorCode TSStep_RosW(TS ts) 1127 { 1128 TS_RosW *ros = (TS_RosW *)ts->data; 1129 RosWTableau tab = ros->tableau; 1130 const PetscInt s = tab->s; 1131 const PetscReal *At = tab->At, *Gamma = tab->Gamma, *ASum = tab->ASum, *GammaInv = tab->GammaInv; 1132 const PetscReal *GammaExplicitCorr = tab->GammaExplicitCorr; 1133 const PetscBool *GammaZeroDiag = tab->GammaZeroDiag; 1134 PetscScalar *w = ros->work; 1135 Vec *Y = ros->Y, Ydot = ros->Ydot, Zdot = ros->Zdot, Zstage = ros->Zstage; 1136 SNES snes; 1137 TSAdapt adapt; 1138 PetscInt i, j, its, lits; 1139 PetscInt rejections = 0; 1140 PetscBool stageok, accept = PETSC_TRUE; 1141 PetscReal next_time_step = ts->time_step; 1142 PetscInt lag; 1143 1144 PetscFunctionBegin; 1145 if (!ts->steprollback) PetscCall(VecCopy(ts->vec_sol, ros->vec_sol_prev)); 1146 1147 ros->status = TS_STEP_INCOMPLETE; 1148 while (!ts->reason && ros->status != TS_STEP_COMPLETE) { 1149 const PetscReal h = ts->time_step; 1150 for (i = 0; i < s; i++) { 1151 ros->stage_time = ts->ptime + h * ASum[i]; 1152 PetscCall(TSPreStage(ts, ros->stage_time)); 1153 if (GammaZeroDiag[i]) { 1154 ros->stage_explicit = PETSC_TRUE; 1155 ros->scoeff = 1.; 1156 } else { 1157 ros->stage_explicit = PETSC_FALSE; 1158 ros->scoeff = 1. / Gamma[i * s + i]; 1159 } 1160 1161 PetscCall(VecCopy(ts->vec_sol, Zstage)); 1162 for (j = 0; j < i; j++) w[j] = At[i * s + j]; 1163 PetscCall(VecMAXPY(Zstage, i, w, Y)); 1164 1165 for (j = 0; j < i; j++) w[j] = 1. / h * GammaInv[i * s + j]; 1166 PetscCall(VecZeroEntries(Zdot)); 1167 PetscCall(VecMAXPY(Zdot, i, w, Y)); 1168 1169 /* Initial guess taken from last stage */ 1170 PetscCall(VecZeroEntries(Y[i])); 1171 1172 if (!ros->stage_explicit) { 1173 PetscCall(TSGetSNES(ts, &snes)); 1174 if (!ros->recompute_jacobian && !i) { 1175 PetscCall(SNESGetLagJacobian(snes, &lag)); 1176 if (lag == 1) { /* use did not set a nontrivial lag, so lag over all stages */ 1177 PetscCall(SNESSetLagJacobian(snes, -2)); /* Recompute the Jacobian on this solve, but not again for the rest of the stages */ 1178 } 1179 } 1180 PetscCall(SNESSolve(snes, NULL, Y[i])); 1181 if (!ros->recompute_jacobian && i == s - 1 && lag == 1) { PetscCall(SNESSetLagJacobian(snes, lag)); /* Set lag back to 1 so we know user did not set it */ } 1182 PetscCall(SNESGetIterationNumber(snes, &its)); 1183 PetscCall(SNESGetLinearSolveIterations(snes, &lits)); 1184 ts->snes_its += its; 1185 ts->ksp_its += lits; 1186 } else { 1187 Mat J, Jp; 1188 PetscCall(VecZeroEntries(Ydot)); /* Evaluate Y[i]=G(t,Ydot=0,Zstage) */ 1189 PetscCall(TSComputeIFunction(ts, ros->stage_time, Zstage, Ydot, Y[i], PETSC_FALSE)); 1190 PetscCall(VecScale(Y[i], -1.0)); 1191 PetscCall(VecAXPY(Y[i], -1.0, Zdot)); /*Y[i] = F(Zstage)-Zdot[=GammaInv*Y]*/ 1192 1193 PetscCall(VecZeroEntries(Zstage)); /* Zstage = GammaExplicitCorr[i,j] * Y[j] */ 1194 for (j = 0; j < i; j++) w[j] = GammaExplicitCorr[i * s + j]; 1195 PetscCall(VecMAXPY(Zstage, i, w, Y)); 1196 1197 /* Y[i] = Y[i] + Jac*Zstage[=Jac*GammaExplicitCorr[i,j] * Y[j]] */ 1198 PetscCall(TSGetIJacobian(ts, &J, &Jp, NULL, NULL)); 1199 PetscCall(TSComputeIJacobian(ts, ros->stage_time, ts->vec_sol, Ydot, 0, J, Jp, PETSC_FALSE)); 1200 PetscCall(MatMult(J, Zstage, Zdot)); 1201 PetscCall(VecAXPY(Y[i], -1.0, Zdot)); 1202 ts->ksp_its += 1; 1203 1204 PetscCall(VecScale(Y[i], h)); 1205 } 1206 PetscCall(TSPostStage(ts, ros->stage_time, i, Y)); 1207 PetscCall(TSGetAdapt(ts, &adapt)); 1208 PetscCall(TSAdaptCheckStage(adapt, ts, ros->stage_time, Y[i], &stageok)); 1209 if (!stageok) goto reject_step; 1210 } 1211 1212 ros->status = TS_STEP_INCOMPLETE; 1213 PetscCall(TSEvaluateStep_RosW(ts, tab->order, ts->vec_sol, NULL)); 1214 ros->status = TS_STEP_PENDING; 1215 PetscCall(TSGetAdapt(ts, &adapt)); 1216 PetscCall(TSAdaptCandidatesClear(adapt)); 1217 PetscCall(TSAdaptCandidateAdd(adapt, tab->name, tab->order, 1, tab->ccfl, (PetscReal)tab->s, PETSC_TRUE)); 1218 PetscCall(TSAdaptChoose(adapt, ts, ts->time_step, NULL, &next_time_step, &accept)); 1219 ros->status = accept ? TS_STEP_COMPLETE : TS_STEP_INCOMPLETE; 1220 if (!accept) { /* Roll back the current step */ 1221 PetscCall(VecCopy(ts->vec_sol0, ts->vec_sol)); 1222 ts->time_step = next_time_step; 1223 goto reject_step; 1224 } 1225 1226 ts->ptime += ts->time_step; 1227 ts->time_step = next_time_step; 1228 break; 1229 1230 reject_step: 1231 ts->reject++; 1232 accept = PETSC_FALSE; 1233 if (!ts->reason && ++rejections > ts->max_reject && ts->max_reject >= 0) { 1234 ts->reason = TS_DIVERGED_STEP_REJECTED; 1235 PetscCall(PetscInfo(ts, "Step=%" PetscInt_FMT ", step rejections %" PetscInt_FMT " greater than current TS allowed, stopping solve\n", ts->steps, rejections)); 1236 } 1237 } 1238 PetscFunctionReturn(PETSC_SUCCESS); 1239 } 1240 1241 static PetscErrorCode TSInterpolate_RosW(TS ts, PetscReal itime, Vec U) 1242 { 1243 TS_RosW *ros = (TS_RosW *)ts->data; 1244 PetscInt s = ros->tableau->s, pinterp = ros->tableau->pinterp, i, j; 1245 PetscReal h; 1246 PetscReal tt, t; 1247 PetscScalar *bt; 1248 const PetscReal *Bt = ros->tableau->binterpt; 1249 const PetscReal *GammaInv = ros->tableau->GammaInv; 1250 PetscScalar *w = ros->work; 1251 Vec *Y = ros->Y; 1252 1253 PetscFunctionBegin; 1254 PetscCheck(Bt, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "TSRosW %s does not have an interpolation formula", ros->tableau->name); 1255 1256 switch (ros->status) { 1257 case TS_STEP_INCOMPLETE: 1258 case TS_STEP_PENDING: 1259 h = ts->time_step; 1260 t = (itime - ts->ptime) / h; 1261 break; 1262 case TS_STEP_COMPLETE: 1263 h = ts->ptime - ts->ptime_prev; 1264 t = (itime - ts->ptime) / h + 1; /* In the interval [0,1] */ 1265 break; 1266 default: 1267 SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_PLIB, "Invalid TSStepStatus"); 1268 } 1269 PetscCall(PetscMalloc1(s, &bt)); 1270 for (i = 0; i < s; i++) bt[i] = 0; 1271 for (j = 0, tt = t; j < pinterp; j++, tt *= t) { 1272 for (i = 0; i < s; i++) bt[i] += Bt[i * pinterp + j] * tt; 1273 } 1274 1275 /* y(t+tt*h) = y(t) + Sum bt(tt) * GammaInv * Ydot */ 1276 /* U <- 0*/ 1277 PetscCall(VecZeroEntries(U)); 1278 /* U <- Sum bt_i * GammaInv(i,1:i) * Y(1:i) */ 1279 for (j = 0; j < s; j++) w[j] = 0; 1280 for (j = 0; j < s; j++) { 1281 for (i = j; i < s; i++) w[j] += bt[i] * GammaInv[i * s + j]; 1282 } 1283 PetscCall(VecMAXPY(U, i, w, Y)); 1284 /* U <- y(t) + U */ 1285 PetscCall(VecAXPY(U, 1, ros->vec_sol_prev)); 1286 1287 PetscCall(PetscFree(bt)); 1288 PetscFunctionReturn(PETSC_SUCCESS); 1289 } 1290 1291 /*------------------------------------------------------------*/ 1292 1293 static PetscErrorCode TSRosWTableauReset(TS ts) 1294 { 1295 TS_RosW *ros = (TS_RosW *)ts->data; 1296 RosWTableau tab = ros->tableau; 1297 1298 PetscFunctionBegin; 1299 if (!tab) PetscFunctionReturn(PETSC_SUCCESS); 1300 PetscCall(VecDestroyVecs(tab->s, &ros->Y)); 1301 PetscCall(PetscFree(ros->work)); 1302 PetscFunctionReturn(PETSC_SUCCESS); 1303 } 1304 1305 static PetscErrorCode TSReset_RosW(TS ts) 1306 { 1307 TS_RosW *ros = (TS_RosW *)ts->data; 1308 1309 PetscFunctionBegin; 1310 PetscCall(TSRosWTableauReset(ts)); 1311 PetscCall(VecDestroy(&ros->Ydot)); 1312 PetscCall(VecDestroy(&ros->Ystage)); 1313 PetscCall(VecDestroy(&ros->Zdot)); 1314 PetscCall(VecDestroy(&ros->Zstage)); 1315 PetscCall(VecDestroy(&ros->vec_sol_prev)); 1316 PetscFunctionReturn(PETSC_SUCCESS); 1317 } 1318 1319 static PetscErrorCode TSRosWGetVecs(TS ts, DM dm, Vec *Ydot, Vec *Zdot, Vec *Ystage, Vec *Zstage) 1320 { 1321 TS_RosW *rw = (TS_RosW *)ts->data; 1322 1323 PetscFunctionBegin; 1324 if (Ydot) { 1325 if (dm && dm != ts->dm) { 1326 PetscCall(DMGetNamedGlobalVector(dm, "TSRosW_Ydot", Ydot)); 1327 } else *Ydot = rw->Ydot; 1328 } 1329 if (Zdot) { 1330 if (dm && dm != ts->dm) { 1331 PetscCall(DMGetNamedGlobalVector(dm, "TSRosW_Zdot", Zdot)); 1332 } else *Zdot = rw->Zdot; 1333 } 1334 if (Ystage) { 1335 if (dm && dm != ts->dm) { 1336 PetscCall(DMGetNamedGlobalVector(dm, "TSRosW_Ystage", Ystage)); 1337 } else *Ystage = rw->Ystage; 1338 } 1339 if (Zstage) { 1340 if (dm && dm != ts->dm) { 1341 PetscCall(DMGetNamedGlobalVector(dm, "TSRosW_Zstage", Zstage)); 1342 } else *Zstage = rw->Zstage; 1343 } 1344 PetscFunctionReturn(PETSC_SUCCESS); 1345 } 1346 1347 static PetscErrorCode TSRosWRestoreVecs(TS ts, DM dm, Vec *Ydot, Vec *Zdot, Vec *Ystage, Vec *Zstage) 1348 { 1349 PetscFunctionBegin; 1350 if (Ydot) { 1351 if (dm && dm != ts->dm) PetscCall(DMRestoreNamedGlobalVector(dm, "TSRosW_Ydot", Ydot)); 1352 } 1353 if (Zdot) { 1354 if (dm && dm != ts->dm) PetscCall(DMRestoreNamedGlobalVector(dm, "TSRosW_Zdot", Zdot)); 1355 } 1356 if (Ystage) { 1357 if (dm && dm != ts->dm) PetscCall(DMRestoreNamedGlobalVector(dm, "TSRosW_Ystage", Ystage)); 1358 } 1359 if (Zstage) { 1360 if (dm && dm != ts->dm) PetscCall(DMRestoreNamedGlobalVector(dm, "TSRosW_Zstage", Zstage)); 1361 } 1362 PetscFunctionReturn(PETSC_SUCCESS); 1363 } 1364 1365 static PetscErrorCode DMCoarsenHook_TSRosW(DM fine, DM coarse, void *ctx) 1366 { 1367 PetscFunctionBegin; 1368 PetscFunctionReturn(PETSC_SUCCESS); 1369 } 1370 1371 static PetscErrorCode DMRestrictHook_TSRosW(DM fine, Mat restrct, Vec rscale, Mat inject, DM coarse, void *ctx) 1372 { 1373 TS ts = (TS)ctx; 1374 Vec Ydot, Zdot, Ystage, Zstage; 1375 Vec Ydotc, Zdotc, Ystagec, Zstagec; 1376 1377 PetscFunctionBegin; 1378 PetscCall(TSRosWGetVecs(ts, fine, &Ydot, &Ystage, &Zdot, &Zstage)); 1379 PetscCall(TSRosWGetVecs(ts, coarse, &Ydotc, &Ystagec, &Zdotc, &Zstagec)); 1380 PetscCall(MatRestrict(restrct, Ydot, Ydotc)); 1381 PetscCall(VecPointwiseMult(Ydotc, rscale, Ydotc)); 1382 PetscCall(MatRestrict(restrct, Ystage, Ystagec)); 1383 PetscCall(VecPointwiseMult(Ystagec, rscale, Ystagec)); 1384 PetscCall(MatRestrict(restrct, Zdot, Zdotc)); 1385 PetscCall(VecPointwiseMult(Zdotc, rscale, Zdotc)); 1386 PetscCall(MatRestrict(restrct, Zstage, Zstagec)); 1387 PetscCall(VecPointwiseMult(Zstagec, rscale, Zstagec)); 1388 PetscCall(TSRosWRestoreVecs(ts, fine, &Ydot, &Ystage, &Zdot, &Zstage)); 1389 PetscCall(TSRosWRestoreVecs(ts, coarse, &Ydotc, &Ystagec, &Zdotc, &Zstagec)); 1390 PetscFunctionReturn(PETSC_SUCCESS); 1391 } 1392 1393 static PetscErrorCode DMSubDomainHook_TSRosW(DM fine, DM coarse, void *ctx) 1394 { 1395 PetscFunctionBegin; 1396 PetscFunctionReturn(PETSC_SUCCESS); 1397 } 1398 1399 static PetscErrorCode DMSubDomainRestrictHook_TSRosW(DM dm, VecScatter gscat, VecScatter lscat, DM subdm, void *ctx) 1400 { 1401 TS ts = (TS)ctx; 1402 Vec Ydot, Zdot, Ystage, Zstage; 1403 Vec Ydots, Zdots, Ystages, Zstages; 1404 1405 PetscFunctionBegin; 1406 PetscCall(TSRosWGetVecs(ts, dm, &Ydot, &Ystage, &Zdot, &Zstage)); 1407 PetscCall(TSRosWGetVecs(ts, subdm, &Ydots, &Ystages, &Zdots, &Zstages)); 1408 1409 PetscCall(VecScatterBegin(gscat, Ydot, Ydots, INSERT_VALUES, SCATTER_FORWARD)); 1410 PetscCall(VecScatterEnd(gscat, Ydot, Ydots, INSERT_VALUES, SCATTER_FORWARD)); 1411 1412 PetscCall(VecScatterBegin(gscat, Ystage, Ystages, INSERT_VALUES, SCATTER_FORWARD)); 1413 PetscCall(VecScatterEnd(gscat, Ystage, Ystages, INSERT_VALUES, SCATTER_FORWARD)); 1414 1415 PetscCall(VecScatterBegin(gscat, Zdot, Zdots, INSERT_VALUES, SCATTER_FORWARD)); 1416 PetscCall(VecScatterEnd(gscat, Zdot, Zdots, INSERT_VALUES, SCATTER_FORWARD)); 1417 1418 PetscCall(VecScatterBegin(gscat, Zstage, Zstages, INSERT_VALUES, SCATTER_FORWARD)); 1419 PetscCall(VecScatterEnd(gscat, Zstage, Zstages, INSERT_VALUES, SCATTER_FORWARD)); 1420 1421 PetscCall(TSRosWRestoreVecs(ts, dm, &Ydot, &Ystage, &Zdot, &Zstage)); 1422 PetscCall(TSRosWRestoreVecs(ts, subdm, &Ydots, &Ystages, &Zdots, &Zstages)); 1423 PetscFunctionReturn(PETSC_SUCCESS); 1424 } 1425 1426 static PetscErrorCode SNESTSFormFunction_RosW(SNES snes, Vec U, Vec F, TS ts) 1427 { 1428 TS_RosW *ros = (TS_RosW *)ts->data; 1429 Vec Ydot, Zdot, Ystage, Zstage; 1430 PetscReal shift = ros->scoeff / ts->time_step; 1431 DM dm, dmsave; 1432 1433 PetscFunctionBegin; 1434 PetscCall(SNESGetDM(snes, &dm)); 1435 PetscCall(TSRosWGetVecs(ts, dm, &Ydot, &Zdot, &Ystage, &Zstage)); 1436 PetscCall(VecWAXPY(Ydot, shift, U, Zdot)); /* Ydot = shift*U + Zdot */ 1437 PetscCall(VecWAXPY(Ystage, 1.0, U, Zstage)); /* Ystage = U + Zstage */ 1438 dmsave = ts->dm; 1439 ts->dm = dm; 1440 PetscCall(TSComputeIFunction(ts, ros->stage_time, Ystage, Ydot, F, PETSC_FALSE)); 1441 ts->dm = dmsave; 1442 PetscCall(TSRosWRestoreVecs(ts, dm, &Ydot, &Zdot, &Ystage, &Zstage)); 1443 PetscFunctionReturn(PETSC_SUCCESS); 1444 } 1445 1446 static PetscErrorCode SNESTSFormJacobian_RosW(SNES snes, Vec U, Mat A, Mat B, TS ts) 1447 { 1448 TS_RosW *ros = (TS_RosW *)ts->data; 1449 Vec Ydot, Zdot, Ystage, Zstage; 1450 PetscReal shift = ros->scoeff / ts->time_step; 1451 DM dm, dmsave; 1452 1453 PetscFunctionBegin; 1454 /* ros->Ydot and ros->Ystage have already been computed in SNESTSFormFunction_RosW (SNES guarantees this) */ 1455 PetscCall(SNESGetDM(snes, &dm)); 1456 PetscCall(TSRosWGetVecs(ts, dm, &Ydot, &Zdot, &Ystage, &Zstage)); 1457 dmsave = ts->dm; 1458 ts->dm = dm; 1459 PetscCall(TSComputeIJacobian(ts, ros->stage_time, Ystage, Ydot, shift, A, B, PETSC_TRUE)); 1460 ts->dm = dmsave; 1461 PetscCall(TSRosWRestoreVecs(ts, dm, &Ydot, &Zdot, &Ystage, &Zstage)); 1462 PetscFunctionReturn(PETSC_SUCCESS); 1463 } 1464 1465 static PetscErrorCode TSRosWTableauSetUp(TS ts) 1466 { 1467 TS_RosW *ros = (TS_RosW *)ts->data; 1468 RosWTableau tab = ros->tableau; 1469 1470 PetscFunctionBegin; 1471 PetscCall(VecDuplicateVecs(ts->vec_sol, tab->s, &ros->Y)); 1472 PetscCall(PetscMalloc1(tab->s, &ros->work)); 1473 PetscFunctionReturn(PETSC_SUCCESS); 1474 } 1475 1476 static PetscErrorCode TSSetUp_RosW(TS ts) 1477 { 1478 TS_RosW *ros = (TS_RosW *)ts->data; 1479 DM dm; 1480 SNES snes; 1481 TSRHSJacobianFn *rhsjacobian; 1482 1483 PetscFunctionBegin; 1484 PetscCall(TSRosWTableauSetUp(ts)); 1485 PetscCall(VecDuplicate(ts->vec_sol, &ros->Ydot)); 1486 PetscCall(VecDuplicate(ts->vec_sol, &ros->Ystage)); 1487 PetscCall(VecDuplicate(ts->vec_sol, &ros->Zdot)); 1488 PetscCall(VecDuplicate(ts->vec_sol, &ros->Zstage)); 1489 PetscCall(VecDuplicate(ts->vec_sol, &ros->vec_sol_prev)); 1490 PetscCall(TSGetDM(ts, &dm)); 1491 PetscCall(DMCoarsenHookAdd(dm, DMCoarsenHook_TSRosW, DMRestrictHook_TSRosW, ts)); 1492 PetscCall(DMSubDomainHookAdd(dm, DMSubDomainHook_TSRosW, DMSubDomainRestrictHook_TSRosW, ts)); 1493 /* Rosenbrock methods are linearly implicit, so set that unless the user has specifically asked for something else */ 1494 PetscCall(TSGetSNES(ts, &snes)); 1495 if (!((PetscObject)snes)->type_name) PetscCall(SNESSetType(snes, SNESKSPONLY)); 1496 PetscCall(DMTSGetRHSJacobian(dm, &rhsjacobian, NULL)); 1497 if (rhsjacobian == TSComputeRHSJacobianConstant) { 1498 Mat Amat, Pmat; 1499 1500 /* Set the SNES matrix to be different from the RHS matrix because there is no way to reconstruct shift*M-J */ 1501 PetscCall(SNESGetJacobian(snes, &Amat, &Pmat, NULL, NULL)); 1502 if (Amat && Amat == ts->Arhs) { 1503 if (Amat == Pmat) { 1504 PetscCall(MatDuplicate(ts->Arhs, MAT_COPY_VALUES, &Amat)); 1505 PetscCall(SNESSetJacobian(snes, Amat, Amat, NULL, NULL)); 1506 } else { 1507 PetscCall(MatDuplicate(ts->Arhs, MAT_COPY_VALUES, &Amat)); 1508 PetscCall(SNESSetJacobian(snes, Amat, NULL, NULL, NULL)); 1509 if (Pmat && Pmat == ts->Brhs) { 1510 PetscCall(MatDuplicate(ts->Brhs, MAT_COPY_VALUES, &Pmat)); 1511 PetscCall(SNESSetJacobian(snes, NULL, Pmat, NULL, NULL)); 1512 PetscCall(MatDestroy(&Pmat)); 1513 } 1514 } 1515 PetscCall(MatDestroy(&Amat)); 1516 } 1517 } 1518 PetscFunctionReturn(PETSC_SUCCESS); 1519 } 1520 /*------------------------------------------------------------*/ 1521 1522 static PetscErrorCode TSSetFromOptions_RosW(TS ts, PetscOptionItems PetscOptionsObject) 1523 { 1524 TS_RosW *ros = (TS_RosW *)ts->data; 1525 SNES snes; 1526 1527 PetscFunctionBegin; 1528 PetscOptionsHeadBegin(PetscOptionsObject, "RosW ODE solver options"); 1529 { 1530 RosWTableauLink link; 1531 PetscInt count, choice; 1532 PetscBool flg; 1533 const char **namelist; 1534 1535 for (link = RosWTableauList, count = 0; link; link = link->next, count++); 1536 PetscCall(PetscMalloc1(count, (char ***)&namelist)); 1537 for (link = RosWTableauList, count = 0; link; link = link->next, count++) namelist[count] = link->tab.name; 1538 PetscCall(PetscOptionsEList("-ts_rosw_type", "Family of Rosenbrock-W method", "TSRosWSetType", (const char *const *)namelist, count, ros->tableau->name, &choice, &flg)); 1539 if (flg) PetscCall(TSRosWSetType(ts, namelist[choice])); 1540 PetscCall(PetscFree(namelist)); 1541 1542 PetscCall(PetscOptionsBool("-ts_rosw_recompute_jacobian", "Recompute the Jacobian at each stage", "TSRosWSetRecomputeJacobian", ros->recompute_jacobian, &ros->recompute_jacobian, NULL)); 1543 } 1544 PetscOptionsHeadEnd(); 1545 /* Rosenbrock methods are linearly implicit, so set that unless the user has specifically asked for something else */ 1546 PetscCall(TSGetSNES(ts, &snes)); 1547 if (!((PetscObject)snes)->type_name) PetscCall(SNESSetType(snes, SNESKSPONLY)); 1548 PetscFunctionReturn(PETSC_SUCCESS); 1549 } 1550 1551 static PetscErrorCode TSView_RosW(TS ts, PetscViewer viewer) 1552 { 1553 TS_RosW *ros = (TS_RosW *)ts->data; 1554 PetscBool iascii; 1555 1556 PetscFunctionBegin; 1557 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii)); 1558 if (iascii) { 1559 RosWTableau tab = ros->tableau; 1560 TSRosWType rostype; 1561 char buf[512]; 1562 PetscInt i; 1563 PetscReal abscissa[512]; 1564 PetscCall(TSRosWGetType(ts, &rostype)); 1565 PetscCall(PetscViewerASCIIPrintf(viewer, " Rosenbrock-W %s\n", rostype)); 1566 PetscCall(PetscFormatRealArray(buf, sizeof(buf), "% 8.6f", tab->s, tab->ASum)); 1567 PetscCall(PetscViewerASCIIPrintf(viewer, " Abscissa of A = %s\n", buf)); 1568 for (i = 0; i < tab->s; i++) abscissa[i] = tab->ASum[i] + tab->Gamma[i]; 1569 PetscCall(PetscFormatRealArray(buf, sizeof(buf), "% 8.6f", tab->s, abscissa)); 1570 PetscCall(PetscViewerASCIIPrintf(viewer, " Abscissa of A+Gamma = %s\n", buf)); 1571 } 1572 PetscFunctionReturn(PETSC_SUCCESS); 1573 } 1574 1575 static PetscErrorCode TSLoad_RosW(TS ts, PetscViewer viewer) 1576 { 1577 SNES snes; 1578 TSAdapt adapt; 1579 1580 PetscFunctionBegin; 1581 PetscCall(TSGetAdapt(ts, &adapt)); 1582 PetscCall(TSAdaptLoad(adapt, viewer)); 1583 PetscCall(TSGetSNES(ts, &snes)); 1584 PetscCall(SNESLoad(snes, viewer)); 1585 /* function and Jacobian context for SNES when used with TS is always ts object */ 1586 PetscCall(SNESSetFunction(snes, NULL, NULL, ts)); 1587 PetscCall(SNESSetJacobian(snes, NULL, NULL, NULL, ts)); 1588 PetscFunctionReturn(PETSC_SUCCESS); 1589 } 1590 1591 /*@ 1592 TSRosWSetType - Set the type of Rosenbrock-W, `TSROSW`, scheme 1593 1594 Logically Collective 1595 1596 Input Parameters: 1597 + ts - timestepping context 1598 - roswtype - type of Rosenbrock-W scheme 1599 1600 Level: beginner 1601 1602 .seealso: [](ch_ts), `TSRosWGetType()`, `TSROSW`, `TSROSW2M`, `TSROSW2P`, `TSROSWRA3PW`, `TSROSWRA34PW2`, `TSROSWRODAS3`, `TSROSWSANDU3`, `TSROSWASSP3P3S1C`, `TSROSWLASSP3P4S2C`, `TSROSWLLSSP3P4S2C`, `TSROSWARK3` 1603 @*/ 1604 PetscErrorCode TSRosWSetType(TS ts, TSRosWType roswtype) 1605 { 1606 PetscFunctionBegin; 1607 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 1608 PetscAssertPointer(roswtype, 2); 1609 PetscTryMethod(ts, "TSRosWSetType_C", (TS, TSRosWType), (ts, roswtype)); 1610 PetscFunctionReturn(PETSC_SUCCESS); 1611 } 1612 1613 /*@ 1614 TSRosWGetType - Get the type of Rosenbrock-W scheme 1615 1616 Logically Collective 1617 1618 Input Parameter: 1619 . ts - timestepping context 1620 1621 Output Parameter: 1622 . rostype - type of Rosenbrock-W scheme 1623 1624 Level: intermediate 1625 1626 .seealso: [](ch_ts), `TSRosWType`, `TSRosWSetType()` 1627 @*/ 1628 PetscErrorCode TSRosWGetType(TS ts, TSRosWType *rostype) 1629 { 1630 PetscFunctionBegin; 1631 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 1632 PetscUseMethod(ts, "TSRosWGetType_C", (TS, TSRosWType *), (ts, rostype)); 1633 PetscFunctionReturn(PETSC_SUCCESS); 1634 } 1635 1636 /*@ 1637 TSRosWSetRecomputeJacobian - Set whether to recompute the Jacobian at each stage. The default is to update the Jacobian once per step. 1638 1639 Logically Collective 1640 1641 Input Parameters: 1642 + ts - timestepping context 1643 - flg - `PETSC_TRUE` to recompute the Jacobian at each stage 1644 1645 Level: intermediate 1646 1647 .seealso: [](ch_ts), `TSRosWType`, `TSRosWGetType()` 1648 @*/ 1649 PetscErrorCode TSRosWSetRecomputeJacobian(TS ts, PetscBool flg) 1650 { 1651 PetscFunctionBegin; 1652 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 1653 PetscTryMethod(ts, "TSRosWSetRecomputeJacobian_C", (TS, PetscBool), (ts, flg)); 1654 PetscFunctionReturn(PETSC_SUCCESS); 1655 } 1656 1657 static PetscErrorCode TSRosWGetType_RosW(TS ts, TSRosWType *rostype) 1658 { 1659 TS_RosW *ros = (TS_RosW *)ts->data; 1660 1661 PetscFunctionBegin; 1662 *rostype = ros->tableau->name; 1663 PetscFunctionReturn(PETSC_SUCCESS); 1664 } 1665 1666 static PetscErrorCode TSRosWSetType_RosW(TS ts, TSRosWType rostype) 1667 { 1668 TS_RosW *ros = (TS_RosW *)ts->data; 1669 PetscBool match; 1670 RosWTableauLink link; 1671 1672 PetscFunctionBegin; 1673 if (ros->tableau) { 1674 PetscCall(PetscStrcmp(ros->tableau->name, rostype, &match)); 1675 if (match) PetscFunctionReturn(PETSC_SUCCESS); 1676 } 1677 for (link = RosWTableauList; link; link = link->next) { 1678 PetscCall(PetscStrcmp(link->tab.name, rostype, &match)); 1679 if (match) { 1680 if (ts->setupcalled) PetscCall(TSRosWTableauReset(ts)); 1681 ros->tableau = &link->tab; 1682 if (ts->setupcalled) PetscCall(TSRosWTableauSetUp(ts)); 1683 ts->default_adapt_type = ros->tableau->bembed ? TSADAPTBASIC : TSADAPTNONE; 1684 PetscFunctionReturn(PETSC_SUCCESS); 1685 } 1686 } 1687 SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_UNKNOWN_TYPE, "Could not find '%s'", rostype); 1688 } 1689 1690 static PetscErrorCode TSRosWSetRecomputeJacobian_RosW(TS ts, PetscBool flg) 1691 { 1692 TS_RosW *ros = (TS_RosW *)ts->data; 1693 1694 PetscFunctionBegin; 1695 ros->recompute_jacobian = flg; 1696 PetscFunctionReturn(PETSC_SUCCESS); 1697 } 1698 1699 static PetscErrorCode TSDestroy_RosW(TS ts) 1700 { 1701 PetscFunctionBegin; 1702 PetscCall(TSReset_RosW(ts)); 1703 if (ts->dm) { 1704 PetscCall(DMCoarsenHookRemove(ts->dm, DMCoarsenHook_TSRosW, DMRestrictHook_TSRosW, ts)); 1705 PetscCall(DMSubDomainHookRemove(ts->dm, DMSubDomainHook_TSRosW, DMSubDomainRestrictHook_TSRosW, ts)); 1706 } 1707 PetscCall(PetscFree(ts->data)); 1708 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSRosWGetType_C", NULL)); 1709 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSRosWSetType_C", NULL)); 1710 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSRosWSetRecomputeJacobian_C", NULL)); 1711 PetscFunctionReturn(PETSC_SUCCESS); 1712 } 1713 1714 /* ------------------------------------------------------------ */ 1715 /*MC 1716 TSROSW - ODE solver using Rosenbrock-W schemes 1717 1718 These methods are intended for problems with well-separated time scales, especially when a slow scale is strongly 1719 nonlinear such that it is expensive to solve with a fully implicit method. The user should provide the stiff part 1720 of the equation using `TSSetIFunction()` and the non-stiff part with `TSSetRHSFunction()`. 1721 1722 Level: beginner 1723 1724 Notes: 1725 This method currently only works with autonomous ODE and DAE. 1726 1727 Consider trying `TSARKIMEX` if the stiff part is strongly nonlinear. 1728 1729 Since this uses a single linear solve per time-step if you wish to lag the jacobian or preconditioner computation you must use also -snes_lag_jacobian_persists true or -snes_lag_jacobian_preconditioner true 1730 1731 Developer Notes: 1732 Rosenbrock-W methods are typically specified for autonomous ODE 1733 .vb 1734 udot = f(u) 1735 .ve 1736 by the stage equations 1737 .vb 1738 k_i = h f(u_0 + sum_j alpha_ij k_j) + h J sum_j gamma_ij k_j 1739 .ve 1740 and step completion formula 1741 .vb 1742 u_1 = u_0 + sum_j b_j k_j 1743 .ve 1744 with step size h and coefficients alpha_ij, gamma_ij, and b_i. Implementing the method in this form would require f(u) 1745 and the Jacobian J to be available, in addition to the shifted matrix I - h gamma_ii J. Following Hairer and Wanner, 1746 we define new variables for the stage equations 1747 .vb 1748 y_i = gamma_ij k_j 1749 .ve 1750 The k_j can be recovered because Gamma is invertible. Let C be the lower triangular part of Gamma^{-1} and define 1751 .vb 1752 A = Alpha Gamma^{-1}, bt^T = b^T Gamma^{-1} 1753 .ve 1754 to rewrite the method as 1755 .vb 1756 [M/(h gamma_ii) - J] y_i = f(u_0 + sum_j a_ij y_j) + M sum_j (c_ij/h) y_j 1757 u_1 = u_0 + sum_j bt_j y_j 1758 .ve 1759 1760 where we have introduced the mass matrix M. Continue by defining 1761 .vb 1762 ydot_i = 1/(h gamma_ii) y_i - sum_j (c_ij/h) y_j 1763 .ve 1764 or, more compactly in tensor notation 1765 .vb 1766 Ydot = 1/h (Gamma^{-1} \otimes I) Y . 1767 .ve 1768 Note that Gamma^{-1} is lower triangular. With this definition of Ydot in terms of known quantities and the current 1769 stage y_i, the stage equations reduce to performing one Newton step (typically with a lagged Jacobian) on the 1770 equation 1771 .vb 1772 g(u_0 + sum_j a_ij y_j + y_i, ydot_i) = 0 1773 .ve 1774 with initial guess y_i = 0. 1775 1776 .seealso: [](ch_ts), `TSCreate()`, `TS`, `TSSetType()`, `TSRosWSetType()`, `TSRosWRegister()`, `TSROSWTHETA1`, `TSROSWTHETA2`, `TSROSW2M`, `TSROSW2P`, `TSROSWRA3PW`, `TSROSWRA34PW2`, `TSROSWRODAS3`, 1777 `TSROSWSANDU3`, `TSROSWASSP3P3S1C`, `TSROSWLASSP3P4S2C`, `TSROSWLLSSP3P4S2C`, `TSROSWGRK4T`, `TSROSWSHAMP4`, `TSROSWVELDD4`, `TSROSW4L`, `TSType` 1778 M*/ 1779 PETSC_EXTERN PetscErrorCode TSCreate_RosW(TS ts) 1780 { 1781 TS_RosW *ros; 1782 1783 PetscFunctionBegin; 1784 PetscCall(TSRosWInitializePackage()); 1785 1786 ts->ops->reset = TSReset_RosW; 1787 ts->ops->destroy = TSDestroy_RosW; 1788 ts->ops->view = TSView_RosW; 1789 ts->ops->load = TSLoad_RosW; 1790 ts->ops->setup = TSSetUp_RosW; 1791 ts->ops->step = TSStep_RosW; 1792 ts->ops->interpolate = TSInterpolate_RosW; 1793 ts->ops->evaluatestep = TSEvaluateStep_RosW; 1794 ts->ops->setfromoptions = TSSetFromOptions_RosW; 1795 ts->ops->snesfunction = SNESTSFormFunction_RosW; 1796 ts->ops->snesjacobian = SNESTSFormJacobian_RosW; 1797 1798 ts->usessnes = PETSC_TRUE; 1799 1800 PetscCall(PetscNew(&ros)); 1801 ts->data = (void *)ros; 1802 1803 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSRosWGetType_C", TSRosWGetType_RosW)); 1804 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSRosWSetType_C", TSRosWSetType_RosW)); 1805 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSRosWSetRecomputeJacobian_C", TSRosWSetRecomputeJacobian_RosW)); 1806 1807 PetscCall(TSRosWSetType(ts, TSRosWDefault)); 1808 PetscFunctionReturn(PETSC_SUCCESS); 1809 } 1810