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