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 PetscCheck(done, 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.", 1120 tab->name, tab->order, order); 1121 *done = PETSC_FALSE; 1122 PetscFunctionReturn(PETSC_SUCCESS); 1123 } 1124 1125 static PetscErrorCode TSStep_RosW(TS ts) 1126 { 1127 TS_RosW *ros = (TS_RosW *)ts->data; 1128 RosWTableau tab = ros->tableau; 1129 const PetscInt s = tab->s; 1130 const PetscReal *At = tab->At, *Gamma = tab->Gamma, *ASum = tab->ASum, *GammaInv = tab->GammaInv; 1131 const PetscReal *GammaExplicitCorr = tab->GammaExplicitCorr; 1132 const PetscBool *GammaZeroDiag = tab->GammaZeroDiag; 1133 PetscScalar *w = ros->work; 1134 Vec *Y = ros->Y, Ydot = ros->Ydot, Zdot = ros->Zdot, Zstage = ros->Zstage; 1135 SNES snes; 1136 TSAdapt adapt; 1137 PetscInt i, j, its, lits; 1138 PetscInt rejections = 0; 1139 PetscBool stageok, accept = PETSC_TRUE; 1140 PetscReal next_time_step = ts->time_step; 1141 PetscInt lag; 1142 1143 PetscFunctionBegin; 1144 if (!ts->steprollback) PetscCall(VecCopy(ts->vec_sol, ros->vec_sol_prev)); 1145 1146 ros->status = TS_STEP_INCOMPLETE; 1147 while (!ts->reason && ros->status != TS_STEP_COMPLETE) { 1148 const PetscReal h = ts->time_step; 1149 for (i = 0; i < s; i++) { 1150 ros->stage_time = ts->ptime + h * ASum[i]; 1151 PetscCall(TSPreStage(ts, ros->stage_time)); 1152 if (GammaZeroDiag[i]) { 1153 ros->stage_explicit = PETSC_TRUE; 1154 ros->scoeff = 1.; 1155 } else { 1156 ros->stage_explicit = PETSC_FALSE; 1157 ros->scoeff = 1. / Gamma[i * s + i]; 1158 } 1159 1160 PetscCall(VecCopy(ts->vec_sol, Zstage)); 1161 for (j = 0; j < i; j++) w[j] = At[i * s + j]; 1162 PetscCall(VecMAXPY(Zstage, i, w, Y)); 1163 1164 for (j = 0; j < i; j++) w[j] = 1. / h * GammaInv[i * s + j]; 1165 PetscCall(VecZeroEntries(Zdot)); 1166 PetscCall(VecMAXPY(Zdot, i, w, Y)); 1167 1168 /* Initial guess taken from last stage */ 1169 PetscCall(VecZeroEntries(Y[i])); 1170 1171 if (!ros->stage_explicit) { 1172 PetscCall(TSGetSNES(ts, &snes)); 1173 if (!ros->recompute_jacobian && !i) { 1174 PetscCall(SNESGetLagJacobian(snes, &lag)); 1175 if (lag == 1) { /* use did not set a nontrivial lag, so lag over all stages */ 1176 PetscCall(SNESSetLagJacobian(snes, -2)); /* Recompute the Jacobian on this solve, but not again for the rest of the stages */ 1177 } 1178 } 1179 PetscCall(SNESSolve(snes, NULL, Y[i])); 1180 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 */ } 1181 PetscCall(SNESGetIterationNumber(snes, &its)); 1182 PetscCall(SNESGetLinearSolveIterations(snes, &lits)); 1183 ts->snes_its += its; 1184 ts->ksp_its += lits; 1185 } else { 1186 Mat J, Jp; 1187 PetscCall(VecZeroEntries(Ydot)); /* Evaluate Y[i]=G(t,Ydot=0,Zstage) */ 1188 PetscCall(TSComputeIFunction(ts, ros->stage_time, Zstage, Ydot, Y[i], PETSC_FALSE)); 1189 PetscCall(VecScale(Y[i], -1.0)); 1190 PetscCall(VecAXPY(Y[i], -1.0, Zdot)); /*Y[i] = F(Zstage)-Zdot[=GammaInv*Y]*/ 1191 1192 PetscCall(VecZeroEntries(Zstage)); /* Zstage = GammaExplicitCorr[i,j] * Y[j] */ 1193 for (j = 0; j < i; j++) w[j] = GammaExplicitCorr[i * s + j]; 1194 PetscCall(VecMAXPY(Zstage, i, w, Y)); 1195 1196 /* Y[i] = Y[i] + Jac*Zstage[=Jac*GammaExplicitCorr[i,j] * Y[j]] */ 1197 PetscCall(TSGetIJacobian(ts, &J, &Jp, NULL, NULL)); 1198 PetscCall(TSComputeIJacobian(ts, ros->stage_time, ts->vec_sol, Ydot, 0, J, Jp, PETSC_FALSE)); 1199 PetscCall(MatMult(J, Zstage, Zdot)); 1200 PetscCall(VecAXPY(Y[i], -1.0, Zdot)); 1201 ts->ksp_its += 1; 1202 1203 PetscCall(VecScale(Y[i], h)); 1204 } 1205 PetscCall(TSPostStage(ts, ros->stage_time, i, Y)); 1206 PetscCall(TSGetAdapt(ts, &adapt)); 1207 PetscCall(TSAdaptCheckStage(adapt, ts, ros->stage_time, Y[i], &stageok)); 1208 if (!stageok) goto reject_step; 1209 } 1210 1211 ros->status = TS_STEP_INCOMPLETE; 1212 PetscCall(TSEvaluateStep_RosW(ts, tab->order, ts->vec_sol, NULL)); 1213 ros->status = TS_STEP_PENDING; 1214 PetscCall(TSGetAdapt(ts, &adapt)); 1215 PetscCall(TSAdaptCandidatesClear(adapt)); 1216 PetscCall(TSAdaptCandidateAdd(adapt, tab->name, tab->order, 1, tab->ccfl, (PetscReal)tab->s, PETSC_TRUE)); 1217 PetscCall(TSAdaptChoose(adapt, ts, ts->time_step, NULL, &next_time_step, &accept)); 1218 ros->status = accept ? TS_STEP_COMPLETE : TS_STEP_INCOMPLETE; 1219 if (!accept) { /* Roll back the current step */ 1220 PetscCall(VecCopy(ts->vec_sol0, ts->vec_sol)); 1221 ts->time_step = next_time_step; 1222 goto reject_step; 1223 } 1224 1225 ts->ptime += ts->time_step; 1226 ts->time_step = next_time_step; 1227 break; 1228 1229 reject_step: 1230 ts->reject++; 1231 accept = PETSC_FALSE; 1232 if (!ts->reason && ++rejections > ts->max_reject && ts->max_reject >= 0) { 1233 ts->reason = TS_DIVERGED_STEP_REJECTED; 1234 PetscCall(PetscInfo(ts, "Step=%" PetscInt_FMT ", step rejections %" PetscInt_FMT " greater than current TS allowed, stopping solve\n", ts->steps, rejections)); 1235 } 1236 } 1237 PetscFunctionReturn(PETSC_SUCCESS); 1238 } 1239 1240 static PetscErrorCode TSInterpolate_RosW(TS ts, PetscReal itime, Vec U) 1241 { 1242 TS_RosW *ros = (TS_RosW *)ts->data; 1243 PetscInt s = ros->tableau->s, pinterp = ros->tableau->pinterp, i, j; 1244 PetscReal h; 1245 PetscReal tt, t; 1246 PetscScalar *bt; 1247 const PetscReal *Bt = ros->tableau->binterpt; 1248 const PetscReal *GammaInv = ros->tableau->GammaInv; 1249 PetscScalar *w = ros->work; 1250 Vec *Y = ros->Y; 1251 1252 PetscFunctionBegin; 1253 PetscCheck(Bt, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "TSRosW %s does not have an interpolation formula", ros->tableau->name); 1254 1255 switch (ros->status) { 1256 case TS_STEP_INCOMPLETE: 1257 case TS_STEP_PENDING: 1258 h = ts->time_step; 1259 t = (itime - ts->ptime) / h; 1260 break; 1261 case TS_STEP_COMPLETE: 1262 h = ts->ptime - ts->ptime_prev; 1263 t = (itime - ts->ptime) / h + 1; /* In the interval [0,1] */ 1264 break; 1265 default: 1266 SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_PLIB, "Invalid TSStepStatus"); 1267 } 1268 PetscCall(PetscMalloc1(s, &bt)); 1269 for (i = 0; i < s; i++) bt[i] = 0; 1270 for (j = 0, tt = t; j < pinterp; j++, tt *= t) { 1271 for (i = 0; i < s; i++) bt[i] += Bt[i * pinterp + j] * tt; 1272 } 1273 1274 /* y(t+tt*h) = y(t) + Sum bt(tt) * GammaInv * Ydot */ 1275 /* U <- 0*/ 1276 PetscCall(VecZeroEntries(U)); 1277 /* U <- Sum bt_i * GammaInv(i,1:i) * Y(1:i) */ 1278 for (j = 0; j < s; j++) w[j] = 0; 1279 for (j = 0; j < s; j++) { 1280 for (i = j; i < s; i++) w[j] += bt[i] * GammaInv[i * s + j]; 1281 } 1282 PetscCall(VecMAXPY(U, i, w, Y)); 1283 /* U <- y(t) + U */ 1284 PetscCall(VecAXPY(U, 1, ros->vec_sol_prev)); 1285 1286 PetscCall(PetscFree(bt)); 1287 PetscFunctionReturn(PETSC_SUCCESS); 1288 } 1289 1290 /*------------------------------------------------------------*/ 1291 1292 static PetscErrorCode TSRosWTableauReset(TS ts) 1293 { 1294 TS_RosW *ros = (TS_RosW *)ts->data; 1295 RosWTableau tab = ros->tableau; 1296 1297 PetscFunctionBegin; 1298 if (!tab) PetscFunctionReturn(PETSC_SUCCESS); 1299 PetscCall(VecDestroyVecs(tab->s, &ros->Y)); 1300 PetscCall(PetscFree(ros->work)); 1301 PetscFunctionReturn(PETSC_SUCCESS); 1302 } 1303 1304 static PetscErrorCode TSReset_RosW(TS ts) 1305 { 1306 TS_RosW *ros = (TS_RosW *)ts->data; 1307 1308 PetscFunctionBegin; 1309 PetscCall(TSRosWTableauReset(ts)); 1310 PetscCall(VecDestroy(&ros->Ydot)); 1311 PetscCall(VecDestroy(&ros->Ystage)); 1312 PetscCall(VecDestroy(&ros->Zdot)); 1313 PetscCall(VecDestroy(&ros->Zstage)); 1314 PetscCall(VecDestroy(&ros->vec_sol_prev)); 1315 PetscFunctionReturn(PETSC_SUCCESS); 1316 } 1317 1318 static PetscErrorCode TSRosWGetVecs(TS ts, DM dm, Vec *Ydot, Vec *Zdot, Vec *Ystage, Vec *Zstage) 1319 { 1320 TS_RosW *rw = (TS_RosW *)ts->data; 1321 1322 PetscFunctionBegin; 1323 if (Ydot) { 1324 if (dm && dm != ts->dm) { 1325 PetscCall(DMGetNamedGlobalVector(dm, "TSRosW_Ydot", Ydot)); 1326 } else *Ydot = rw->Ydot; 1327 } 1328 if (Zdot) { 1329 if (dm && dm != ts->dm) { 1330 PetscCall(DMGetNamedGlobalVector(dm, "TSRosW_Zdot", Zdot)); 1331 } else *Zdot = rw->Zdot; 1332 } 1333 if (Ystage) { 1334 if (dm && dm != ts->dm) { 1335 PetscCall(DMGetNamedGlobalVector(dm, "TSRosW_Ystage", Ystage)); 1336 } else *Ystage = rw->Ystage; 1337 } 1338 if (Zstage) { 1339 if (dm && dm != ts->dm) { 1340 PetscCall(DMGetNamedGlobalVector(dm, "TSRosW_Zstage", Zstage)); 1341 } else *Zstage = rw->Zstage; 1342 } 1343 PetscFunctionReturn(PETSC_SUCCESS); 1344 } 1345 1346 static PetscErrorCode TSRosWRestoreVecs(TS ts, DM dm, Vec *Ydot, Vec *Zdot, Vec *Ystage, Vec *Zstage) 1347 { 1348 PetscFunctionBegin; 1349 if (Ydot) { 1350 if (dm && dm != ts->dm) PetscCall(DMRestoreNamedGlobalVector(dm, "TSRosW_Ydot", Ydot)); 1351 } 1352 if (Zdot) { 1353 if (dm && dm != ts->dm) PetscCall(DMRestoreNamedGlobalVector(dm, "TSRosW_Zdot", Zdot)); 1354 } 1355 if (Ystage) { 1356 if (dm && dm != ts->dm) PetscCall(DMRestoreNamedGlobalVector(dm, "TSRosW_Ystage", Ystage)); 1357 } 1358 if (Zstage) { 1359 if (dm && dm != ts->dm) PetscCall(DMRestoreNamedGlobalVector(dm, "TSRosW_Zstage", Zstage)); 1360 } 1361 PetscFunctionReturn(PETSC_SUCCESS); 1362 } 1363 1364 static PetscErrorCode DMCoarsenHook_TSRosW(DM fine, DM coarse, void *ctx) 1365 { 1366 PetscFunctionBegin; 1367 PetscFunctionReturn(PETSC_SUCCESS); 1368 } 1369 1370 static PetscErrorCode DMRestrictHook_TSRosW(DM fine, Mat restrct, Vec rscale, Mat inject, DM coarse, void *ctx) 1371 { 1372 TS ts = (TS)ctx; 1373 Vec Ydot, Zdot, Ystage, Zstage; 1374 Vec Ydotc, Zdotc, Ystagec, Zstagec; 1375 1376 PetscFunctionBegin; 1377 PetscCall(TSRosWGetVecs(ts, fine, &Ydot, &Ystage, &Zdot, &Zstage)); 1378 PetscCall(TSRosWGetVecs(ts, coarse, &Ydotc, &Ystagec, &Zdotc, &Zstagec)); 1379 PetscCall(MatRestrict(restrct, Ydot, Ydotc)); 1380 PetscCall(VecPointwiseMult(Ydotc, rscale, Ydotc)); 1381 PetscCall(MatRestrict(restrct, Ystage, Ystagec)); 1382 PetscCall(VecPointwiseMult(Ystagec, rscale, Ystagec)); 1383 PetscCall(MatRestrict(restrct, Zdot, Zdotc)); 1384 PetscCall(VecPointwiseMult(Zdotc, rscale, Zdotc)); 1385 PetscCall(MatRestrict(restrct, Zstage, Zstagec)); 1386 PetscCall(VecPointwiseMult(Zstagec, rscale, Zstagec)); 1387 PetscCall(TSRosWRestoreVecs(ts, fine, &Ydot, &Ystage, &Zdot, &Zstage)); 1388 PetscCall(TSRosWRestoreVecs(ts, coarse, &Ydotc, &Ystagec, &Zdotc, &Zstagec)); 1389 PetscFunctionReturn(PETSC_SUCCESS); 1390 } 1391 1392 static PetscErrorCode DMSubDomainHook_TSRosW(DM fine, DM coarse, void *ctx) 1393 { 1394 PetscFunctionBegin; 1395 PetscFunctionReturn(PETSC_SUCCESS); 1396 } 1397 1398 static PetscErrorCode DMSubDomainRestrictHook_TSRosW(DM dm, VecScatter gscat, VecScatter lscat, DM subdm, void *ctx) 1399 { 1400 TS ts = (TS)ctx; 1401 Vec Ydot, Zdot, Ystage, Zstage; 1402 Vec Ydots, Zdots, Ystages, Zstages; 1403 1404 PetscFunctionBegin; 1405 PetscCall(TSRosWGetVecs(ts, dm, &Ydot, &Ystage, &Zdot, &Zstage)); 1406 PetscCall(TSRosWGetVecs(ts, subdm, &Ydots, &Ystages, &Zdots, &Zstages)); 1407 1408 PetscCall(VecScatterBegin(gscat, Ydot, Ydots, INSERT_VALUES, SCATTER_FORWARD)); 1409 PetscCall(VecScatterEnd(gscat, Ydot, Ydots, INSERT_VALUES, SCATTER_FORWARD)); 1410 1411 PetscCall(VecScatterBegin(gscat, Ystage, Ystages, INSERT_VALUES, SCATTER_FORWARD)); 1412 PetscCall(VecScatterEnd(gscat, Ystage, Ystages, INSERT_VALUES, SCATTER_FORWARD)); 1413 1414 PetscCall(VecScatterBegin(gscat, Zdot, Zdots, INSERT_VALUES, SCATTER_FORWARD)); 1415 PetscCall(VecScatterEnd(gscat, Zdot, Zdots, INSERT_VALUES, SCATTER_FORWARD)); 1416 1417 PetscCall(VecScatterBegin(gscat, Zstage, Zstages, INSERT_VALUES, SCATTER_FORWARD)); 1418 PetscCall(VecScatterEnd(gscat, Zstage, Zstages, INSERT_VALUES, SCATTER_FORWARD)); 1419 1420 PetscCall(TSRosWRestoreVecs(ts, dm, &Ydot, &Ystage, &Zdot, &Zstage)); 1421 PetscCall(TSRosWRestoreVecs(ts, subdm, &Ydots, &Ystages, &Zdots, &Zstages)); 1422 PetscFunctionReturn(PETSC_SUCCESS); 1423 } 1424 1425 static PetscErrorCode SNESTSFormFunction_RosW(SNES snes, Vec U, Vec F, TS ts) 1426 { 1427 TS_RosW *ros = (TS_RosW *)ts->data; 1428 Vec Ydot, Zdot, Ystage, Zstage; 1429 PetscReal shift = ros->scoeff / ts->time_step; 1430 DM dm, dmsave; 1431 1432 PetscFunctionBegin; 1433 PetscCall(SNESGetDM(snes, &dm)); 1434 PetscCall(TSRosWGetVecs(ts, dm, &Ydot, &Zdot, &Ystage, &Zstage)); 1435 PetscCall(VecWAXPY(Ydot, shift, U, Zdot)); /* Ydot = shift*U + Zdot */ 1436 PetscCall(VecWAXPY(Ystage, 1.0, U, Zstage)); /* Ystage = U + Zstage */ 1437 dmsave = ts->dm; 1438 ts->dm = dm; 1439 PetscCall(TSComputeIFunction(ts, ros->stage_time, Ystage, Ydot, F, PETSC_FALSE)); 1440 ts->dm = dmsave; 1441 PetscCall(TSRosWRestoreVecs(ts, dm, &Ydot, &Zdot, &Ystage, &Zstage)); 1442 PetscFunctionReturn(PETSC_SUCCESS); 1443 } 1444 1445 static PetscErrorCode SNESTSFormJacobian_RosW(SNES snes, Vec U, Mat A, Mat B, TS ts) 1446 { 1447 TS_RosW *ros = (TS_RosW *)ts->data; 1448 Vec Ydot, Zdot, Ystage, Zstage; 1449 PetscReal shift = ros->scoeff / ts->time_step; 1450 DM dm, dmsave; 1451 1452 PetscFunctionBegin; 1453 /* ros->Ydot and ros->Ystage have already been computed in SNESTSFormFunction_RosW (SNES guarantees this) */ 1454 PetscCall(SNESGetDM(snes, &dm)); 1455 PetscCall(TSRosWGetVecs(ts, dm, &Ydot, &Zdot, &Ystage, &Zstage)); 1456 dmsave = ts->dm; 1457 ts->dm = dm; 1458 PetscCall(TSComputeIJacobian(ts, ros->stage_time, Ystage, Ydot, shift, A, B, PETSC_TRUE)); 1459 ts->dm = dmsave; 1460 PetscCall(TSRosWRestoreVecs(ts, dm, &Ydot, &Zdot, &Ystage, &Zstage)); 1461 PetscFunctionReturn(PETSC_SUCCESS); 1462 } 1463 1464 static PetscErrorCode TSRosWTableauSetUp(TS ts) 1465 { 1466 TS_RosW *ros = (TS_RosW *)ts->data; 1467 RosWTableau tab = ros->tableau; 1468 1469 PetscFunctionBegin; 1470 PetscCall(VecDuplicateVecs(ts->vec_sol, tab->s, &ros->Y)); 1471 PetscCall(PetscMalloc1(tab->s, &ros->work)); 1472 PetscFunctionReturn(PETSC_SUCCESS); 1473 } 1474 1475 static PetscErrorCode TSSetUp_RosW(TS ts) 1476 { 1477 TS_RosW *ros = (TS_RosW *)ts->data; 1478 DM dm; 1479 SNES snes; 1480 TSRHSJacobianFn *rhsjacobian; 1481 1482 PetscFunctionBegin; 1483 PetscCall(TSRosWTableauSetUp(ts)); 1484 PetscCall(VecDuplicate(ts->vec_sol, &ros->Ydot)); 1485 PetscCall(VecDuplicate(ts->vec_sol, &ros->Ystage)); 1486 PetscCall(VecDuplicate(ts->vec_sol, &ros->Zdot)); 1487 PetscCall(VecDuplicate(ts->vec_sol, &ros->Zstage)); 1488 PetscCall(VecDuplicate(ts->vec_sol, &ros->vec_sol_prev)); 1489 PetscCall(TSGetDM(ts, &dm)); 1490 PetscCall(DMCoarsenHookAdd(dm, DMCoarsenHook_TSRosW, DMRestrictHook_TSRosW, ts)); 1491 PetscCall(DMSubDomainHookAdd(dm, DMSubDomainHook_TSRosW, DMSubDomainRestrictHook_TSRosW, ts)); 1492 /* Rosenbrock methods are linearly implicit, so set that unless the user has specifically asked for something else */ 1493 PetscCall(TSGetSNES(ts, &snes)); 1494 if (!((PetscObject)snes)->type_name) PetscCall(SNESSetType(snes, SNESKSPONLY)); 1495 PetscCall(DMTSGetRHSJacobian(dm, &rhsjacobian, NULL)); 1496 if (rhsjacobian == TSComputeRHSJacobianConstant) { 1497 Mat Amat, Pmat; 1498 1499 /* Set the SNES matrix to be different from the RHS matrix because there is no way to reconstruct shift*M-J */ 1500 PetscCall(SNESGetJacobian(snes, &Amat, &Pmat, NULL, NULL)); 1501 if (Amat && Amat == ts->Arhs) { 1502 if (Amat == Pmat) { 1503 PetscCall(MatDuplicate(ts->Arhs, MAT_COPY_VALUES, &Amat)); 1504 PetscCall(SNESSetJacobian(snes, Amat, Amat, NULL, NULL)); 1505 } else { 1506 PetscCall(MatDuplicate(ts->Arhs, MAT_COPY_VALUES, &Amat)); 1507 PetscCall(SNESSetJacobian(snes, Amat, NULL, NULL, NULL)); 1508 if (Pmat && Pmat == ts->Brhs) { 1509 PetscCall(MatDuplicate(ts->Brhs, MAT_COPY_VALUES, &Pmat)); 1510 PetscCall(SNESSetJacobian(snes, NULL, Pmat, NULL, NULL)); 1511 PetscCall(MatDestroy(&Pmat)); 1512 } 1513 } 1514 PetscCall(MatDestroy(&Amat)); 1515 } 1516 } 1517 PetscFunctionReturn(PETSC_SUCCESS); 1518 } 1519 /*------------------------------------------------------------*/ 1520 1521 static PetscErrorCode TSSetFromOptions_RosW(TS ts, PetscOptionItems PetscOptionsObject) 1522 { 1523 TS_RosW *ros = (TS_RosW *)ts->data; 1524 SNES snes; 1525 1526 PetscFunctionBegin; 1527 PetscOptionsHeadBegin(PetscOptionsObject, "RosW ODE solver options"); 1528 { 1529 RosWTableauLink link; 1530 PetscInt count, choice; 1531 PetscBool flg; 1532 const char **namelist; 1533 1534 for (link = RosWTableauList, count = 0; link; link = link->next, count++); 1535 PetscCall(PetscMalloc1(count, (char ***)&namelist)); 1536 for (link = RosWTableauList, count = 0; link; link = link->next, count++) namelist[count] = link->tab.name; 1537 PetscCall(PetscOptionsEList("-ts_rosw_type", "Family of Rosenbrock-W method", "TSRosWSetType", (const char *const *)namelist, count, ros->tableau->name, &choice, &flg)); 1538 if (flg) PetscCall(TSRosWSetType(ts, namelist[choice])); 1539 PetscCall(PetscFree(namelist)); 1540 1541 PetscCall(PetscOptionsBool("-ts_rosw_recompute_jacobian", "Recompute the Jacobian at each stage", "TSRosWSetRecomputeJacobian", ros->recompute_jacobian, &ros->recompute_jacobian, NULL)); 1542 } 1543 PetscOptionsHeadEnd(); 1544 /* Rosenbrock methods are linearly implicit, so set that unless the user has specifically asked for something else */ 1545 PetscCall(TSGetSNES(ts, &snes)); 1546 if (!((PetscObject)snes)->type_name) PetscCall(SNESSetType(snes, SNESKSPONLY)); 1547 PetscFunctionReturn(PETSC_SUCCESS); 1548 } 1549 1550 static PetscErrorCode TSView_RosW(TS ts, PetscViewer viewer) 1551 { 1552 TS_RosW *ros = (TS_RosW *)ts->data; 1553 PetscBool isascii; 1554 1555 PetscFunctionBegin; 1556 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii)); 1557 if (isascii) { 1558 RosWTableau tab = ros->tableau; 1559 TSRosWType rostype; 1560 char buf[512]; 1561 PetscInt i; 1562 PetscReal abscissa[512]; 1563 PetscCall(TSRosWGetType(ts, &rostype)); 1564 PetscCall(PetscViewerASCIIPrintf(viewer, " Rosenbrock-W %s\n", rostype)); 1565 PetscCall(PetscFormatRealArray(buf, sizeof(buf), "% 8.6f", tab->s, tab->ASum)); 1566 PetscCall(PetscViewerASCIIPrintf(viewer, " Abscissa of A = %s\n", buf)); 1567 for (i = 0; i < tab->s; i++) abscissa[i] = tab->ASum[i] + tab->GammaSum[i]; 1568 PetscCall(PetscFormatRealArray(buf, sizeof(buf), "% 8.6f", tab->s, abscissa)); 1569 PetscCall(PetscViewerASCIIPrintf(viewer, " Abscissa of A+Gamma = %s\n", buf)); 1570 } 1571 PetscFunctionReturn(PETSC_SUCCESS); 1572 } 1573 1574 static PetscErrorCode TSLoad_RosW(TS ts, PetscViewer viewer) 1575 { 1576 SNES snes; 1577 TSAdapt adapt; 1578 1579 PetscFunctionBegin; 1580 PetscCall(TSGetAdapt(ts, &adapt)); 1581 PetscCall(TSAdaptLoad(adapt, viewer)); 1582 PetscCall(TSGetSNES(ts, &snes)); 1583 PetscCall(SNESLoad(snes, viewer)); 1584 /* function and Jacobian context for SNES when used with TS is always ts object */ 1585 PetscCall(SNESSetFunction(snes, NULL, NULL, ts)); 1586 PetscCall(SNESSetJacobian(snes, NULL, NULL, NULL, ts)); 1587 PetscFunctionReturn(PETSC_SUCCESS); 1588 } 1589 1590 /*@ 1591 TSRosWSetType - Set the type of Rosenbrock-W, `TSROSW`, scheme 1592 1593 Logically Collective 1594 1595 Input Parameters: 1596 + ts - timestepping context 1597 - roswtype - type of Rosenbrock-W scheme 1598 1599 Level: beginner 1600 1601 .seealso: [](ch_ts), `TSRosWGetType()`, `TSROSW`, `TSROSW2M`, `TSROSW2P`, `TSROSWRA3PW`, `TSROSWRA34PW2`, `TSROSWRODAS3`, `TSROSWSANDU3`, `TSROSWASSP3P3S1C`, `TSROSWLASSP3P4S2C`, `TSROSWLLSSP3P4S2C`, `TSROSWARK3` 1602 @*/ 1603 PetscErrorCode TSRosWSetType(TS ts, TSRosWType roswtype) 1604 { 1605 PetscFunctionBegin; 1606 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 1607 PetscAssertPointer(roswtype, 2); 1608 PetscTryMethod(ts, "TSRosWSetType_C", (TS, TSRosWType), (ts, roswtype)); 1609 PetscFunctionReturn(PETSC_SUCCESS); 1610 } 1611 1612 /*@ 1613 TSRosWGetType - Get the type of Rosenbrock-W scheme 1614 1615 Logically Collective 1616 1617 Input Parameter: 1618 . ts - timestepping context 1619 1620 Output Parameter: 1621 . rostype - type of Rosenbrock-W scheme 1622 1623 Level: intermediate 1624 1625 .seealso: [](ch_ts), `TSRosWType`, `TSRosWSetType()` 1626 @*/ 1627 PetscErrorCode TSRosWGetType(TS ts, TSRosWType *rostype) 1628 { 1629 PetscFunctionBegin; 1630 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 1631 PetscUseMethod(ts, "TSRosWGetType_C", (TS, TSRosWType *), (ts, rostype)); 1632 PetscFunctionReturn(PETSC_SUCCESS); 1633 } 1634 1635 /*@ 1636 TSRosWSetRecomputeJacobian - Set whether to recompute the Jacobian at each stage. The default is to update the Jacobian once per step. 1637 1638 Logically Collective 1639 1640 Input Parameters: 1641 + ts - timestepping context 1642 - flg - `PETSC_TRUE` to recompute the Jacobian at each stage 1643 1644 Level: intermediate 1645 1646 .seealso: [](ch_ts), `TSRosWType`, `TSRosWGetType()` 1647 @*/ 1648 PetscErrorCode TSRosWSetRecomputeJacobian(TS ts, PetscBool flg) 1649 { 1650 PetscFunctionBegin; 1651 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 1652 PetscTryMethod(ts, "TSRosWSetRecomputeJacobian_C", (TS, PetscBool), (ts, flg)); 1653 PetscFunctionReturn(PETSC_SUCCESS); 1654 } 1655 1656 static PetscErrorCode TSRosWGetType_RosW(TS ts, TSRosWType *rostype) 1657 { 1658 TS_RosW *ros = (TS_RosW *)ts->data; 1659 1660 PetscFunctionBegin; 1661 *rostype = ros->tableau->name; 1662 PetscFunctionReturn(PETSC_SUCCESS); 1663 } 1664 1665 static PetscErrorCode TSRosWSetType_RosW(TS ts, TSRosWType rostype) 1666 { 1667 TS_RosW *ros = (TS_RosW *)ts->data; 1668 PetscBool match; 1669 RosWTableauLink link; 1670 1671 PetscFunctionBegin; 1672 if (ros->tableau) { 1673 PetscCall(PetscStrcmp(ros->tableau->name, rostype, &match)); 1674 if (match) PetscFunctionReturn(PETSC_SUCCESS); 1675 } 1676 for (link = RosWTableauList; link; link = link->next) { 1677 PetscCall(PetscStrcmp(link->tab.name, rostype, &match)); 1678 if (match) { 1679 if (ts->setupcalled) PetscCall(TSRosWTableauReset(ts)); 1680 ros->tableau = &link->tab; 1681 if (ts->setupcalled) PetscCall(TSRosWTableauSetUp(ts)); 1682 ts->default_adapt_type = ros->tableau->bembed ? TSADAPTBASIC : TSADAPTNONE; 1683 PetscFunctionReturn(PETSC_SUCCESS); 1684 } 1685 } 1686 SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_UNKNOWN_TYPE, "Could not find '%s'", rostype); 1687 } 1688 1689 static PetscErrorCode TSRosWSetRecomputeJacobian_RosW(TS ts, PetscBool flg) 1690 { 1691 TS_RosW *ros = (TS_RosW *)ts->data; 1692 1693 PetscFunctionBegin; 1694 ros->recompute_jacobian = flg; 1695 PetscFunctionReturn(PETSC_SUCCESS); 1696 } 1697 1698 static PetscErrorCode TSDestroy_RosW(TS ts) 1699 { 1700 PetscFunctionBegin; 1701 PetscCall(TSReset_RosW(ts)); 1702 if (ts->dm) { 1703 PetscCall(DMCoarsenHookRemove(ts->dm, DMCoarsenHook_TSRosW, DMRestrictHook_TSRosW, ts)); 1704 PetscCall(DMSubDomainHookRemove(ts->dm, DMSubDomainHook_TSRosW, DMSubDomainRestrictHook_TSRosW, ts)); 1705 } 1706 PetscCall(PetscFree(ts->data)); 1707 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSRosWGetType_C", NULL)); 1708 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSRosWSetType_C", NULL)); 1709 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSRosWSetRecomputeJacobian_C", NULL)); 1710 PetscFunctionReturn(PETSC_SUCCESS); 1711 } 1712 1713 /* ------------------------------------------------------------ */ 1714 /*MC 1715 TSROSW - ODE solver using Rosenbrock-W schemes 1716 1717 These methods are intended for problems with well-separated time scales, especially when a slow scale is strongly 1718 nonlinear such that it is expensive to solve with a fully implicit method. The user should provide the stiff part 1719 of the equation using `TSSetIFunction()` and the non-stiff part with `TSSetRHSFunction()`. 1720 1721 Level: beginner 1722 1723 Notes: 1724 This method currently only works with autonomous ODE and DAE. 1725 1726 Consider trying `TSARKIMEX` if the stiff part is strongly nonlinear. 1727 1728 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 1729 1730 Developer Notes: 1731 Rosenbrock-W methods are typically specified for autonomous ODE 1732 .vb 1733 udot = f(u) 1734 .ve 1735 by the stage equations 1736 .vb 1737 k_i = h f(u_0 + sum_j alpha_ij k_j) + h J sum_j gamma_ij k_j 1738 .ve 1739 and step completion formula 1740 .vb 1741 u_1 = u_0 + sum_j b_j k_j 1742 .ve 1743 with step size h and coefficients alpha_ij, gamma_ij, and b_i. Implementing the method in this form would require f(u) 1744 and the Jacobian J to be available, in addition to the shifted matrix I - h gamma_ii J. Following Hairer and Wanner, 1745 we define new variables for the stage equations 1746 .vb 1747 y_i = gamma_ij k_j 1748 .ve 1749 The k_j can be recovered because Gamma is invertible. Let C be the lower triangular part of Gamma^{-1} and define 1750 .vb 1751 A = Alpha Gamma^{-1}, bt^T = b^T Gamma^{-1} 1752 .ve 1753 to rewrite the method as 1754 .vb 1755 [M/(h gamma_ii) - J] y_i = f(u_0 + sum_j a_ij y_j) + M sum_j (c_ij/h) y_j 1756 u_1 = u_0 + sum_j bt_j y_j 1757 .ve 1758 1759 where we have introduced the mass matrix M. Continue by defining 1760 .vb 1761 ydot_i = 1/(h gamma_ii) y_i - sum_j (c_ij/h) y_j 1762 .ve 1763 or, more compactly in tensor notation 1764 .vb 1765 Ydot = 1/h (Gamma^{-1} \otimes I) Y . 1766 .ve 1767 Note that Gamma^{-1} is lower triangular. With this definition of Ydot in terms of known quantities and the current 1768 stage y_i, the stage equations reduce to performing one Newton step (typically with a lagged Jacobian) on the 1769 equation 1770 .vb 1771 g(u_0 + sum_j a_ij y_j + y_i, ydot_i) = 0 1772 .ve 1773 with initial guess y_i = 0. 1774 1775 .seealso: [](ch_ts), `TSCreate()`, `TS`, `TSSetType()`, `TSRosWSetType()`, `TSRosWRegister()`, `TSROSWTHETA1`, `TSROSWTHETA2`, `TSROSW2M`, `TSROSW2P`, `TSROSWRA3PW`, `TSROSWRA34PW2`, `TSROSWRODAS3`, 1776 `TSROSWSANDU3`, `TSROSWASSP3P3S1C`, `TSROSWLASSP3P4S2C`, `TSROSWLLSSP3P4S2C`, `TSROSWGRK4T`, `TSROSWSHAMP4`, `TSROSWVELDD4`, `TSROSW4L`, `TSType` 1777 M*/ 1778 PETSC_EXTERN PetscErrorCode TSCreate_RosW(TS ts) 1779 { 1780 TS_RosW *ros; 1781 1782 PetscFunctionBegin; 1783 PetscCall(TSRosWInitializePackage()); 1784 1785 ts->ops->reset = TSReset_RosW; 1786 ts->ops->destroy = TSDestroy_RosW; 1787 ts->ops->view = TSView_RosW; 1788 ts->ops->load = TSLoad_RosW; 1789 ts->ops->setup = TSSetUp_RosW; 1790 ts->ops->step = TSStep_RosW; 1791 ts->ops->interpolate = TSInterpolate_RosW; 1792 ts->ops->evaluatestep = TSEvaluateStep_RosW; 1793 ts->ops->setfromoptions = TSSetFromOptions_RosW; 1794 ts->ops->snesfunction = SNESTSFormFunction_RosW; 1795 ts->ops->snesjacobian = SNESTSFormJacobian_RosW; 1796 1797 ts->usessnes = PETSC_TRUE; 1798 1799 PetscCall(PetscNew(&ros)); 1800 ts->data = (void *)ros; 1801 1802 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSRosWGetType_C", TSRosWGetType_RosW)); 1803 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSRosWSetType_C", TSRosWSetType_RosW)); 1804 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSRosWSetRecomputeJacobian_C", TSRosWSetRecomputeJacobian_RosW)); 1805 1806 PetscCall(TSRosWSetType(ts, TSRosWDefault)); 1807 PetscFunctionReturn(PETSC_SUCCESS); 1808 } 1809