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