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 MatStructure str; 989 990 PetscFunctionBegin; 991 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 992 next_time_step = ts->time_step; 993 accept = PETSC_TRUE; 994 ros->status = TS_STEP_INCOMPLETE; 995 996 for (reject=0; reject<ts->max_reject && !ts->reason; reject++,ts->reject++) { 997 const PetscReal h = ts->time_step; 998 ierr = TSPreStep(ts);CHKERRQ(ierr); 999 ierr = VecCopy(ts->vec_sol,ros->VecSolPrev);CHKERRQ(ierr); /*move this at the end*/ 1000 for (i=0; i<s; i++) { 1001 ros->stage_time = ts->ptime + h*ASum[i]; 1002 ierr = TSPreStage(ts,ros->stage_time);CHKERRQ(ierr); 1003 if (GammaZeroDiag[i]) { 1004 ros->stage_explicit = PETSC_TRUE; 1005 ros->scoeff = 1.; 1006 } else { 1007 ros->stage_explicit = PETSC_FALSE; 1008 ros->scoeff = 1./Gamma[i*s+i]; 1009 } 1010 1011 ierr = VecCopy(ts->vec_sol,Zstage);CHKERRQ(ierr); 1012 for (j=0; j<i; j++) w[j] = At[i*s+j]; 1013 ierr = VecMAXPY(Zstage,i,w,Y);CHKERRQ(ierr); 1014 1015 for (j=0; j<i; j++) w[j] = 1./h * GammaInv[i*s+j]; 1016 ierr = VecZeroEntries(Zdot);CHKERRQ(ierr); 1017 ierr = VecMAXPY(Zdot,i,w,Y);CHKERRQ(ierr); 1018 1019 /* Initial guess taken from last stage */ 1020 ierr = VecZeroEntries(Y[i]);CHKERRQ(ierr); 1021 1022 if (!ros->stage_explicit) { 1023 if (!ros->recompute_jacobian && !i) { 1024 ierr = SNESSetLagJacobian(snes,-2);CHKERRQ(ierr); /* Recompute the Jacobian on this solve, but not again */ 1025 } 1026 ierr = SNESSolve(snes,NULL,Y[i]);CHKERRQ(ierr); 1027 ierr = SNESGetIterationNumber(snes,&its);CHKERRQ(ierr); 1028 ierr = SNESGetLinearSolveIterations(snes,&lits);CHKERRQ(ierr); 1029 ts->snes_its += its; ts->ksp_its += lits; 1030 ierr = TSGetAdapt(ts,&adapt);CHKERRQ(ierr); 1031 ierr = TSAdaptCheckStage(adapt,ts,&accept);CHKERRQ(ierr); 1032 if (!accept) goto reject_step; 1033 } else { 1034 Mat J,Jp; 1035 ierr = VecZeroEntries(Ydot);CHKERRQ(ierr); /* Evaluate Y[i]=G(t,Ydot=0,Zstage) */ 1036 ierr = TSComputeIFunction(ts,ros->stage_time,Zstage,Ydot,Y[i],PETSC_FALSE);CHKERRQ(ierr); 1037 ierr = VecScale(Y[i],-1.0);CHKERRQ(ierr); 1038 ierr = VecAXPY(Y[i],-1.0,Zdot);CHKERRQ(ierr); /*Y[i]=F(Zstage)-Zdot[=GammaInv*Y]*/ 1039 1040 ierr = VecZeroEntries(Zstage);CHKERRQ(ierr); /* Zstage = GammaExplicitCorr[i,j] * Y[j] */ 1041 for (j=0; j<i; j++) w[j] = GammaExplicitCorr[i*s+j]; 1042 ierr = VecMAXPY(Zstage,i,w,Y);CHKERRQ(ierr); 1043 /*Y[i] += Y[i] + Jac*Zstage[=Jac*GammaExplicitCorr[i,j] * Y[j]] */ 1044 str = SAME_NONZERO_PATTERN; 1045 ierr = TSGetIJacobian(ts,&J,&Jp,NULL,NULL);CHKERRQ(ierr); 1046 ierr = TSComputeIJacobian(ts,ros->stage_time,ts->vec_sol,Ydot,0,&J,&Jp,&str,PETSC_FALSE);CHKERRQ(ierr); 1047 ierr = MatMult(J,Zstage,Zdot);CHKERRQ(ierr); 1048 1049 ierr = VecAXPY(Y[i],-1.0,Zdot);CHKERRQ(ierr); 1050 ierr = VecScale(Y[i],h); 1051 ts->ksp_its += 1; 1052 } 1053 ierr = TSPostStage(ts,ros->stage_time,i,Y);CHKERRQ(ierr); 1054 } 1055 ierr = TSEvaluateStep(ts,tab->order,ts->vec_sol,NULL);CHKERRQ(ierr); 1056 ros->status = TS_STEP_PENDING; 1057 1058 /* Register only the current method as a candidate because we're not supporting multiple candidates yet. */ 1059 ierr = TSGetAdapt(ts,&adapt);CHKERRQ(ierr); 1060 ierr = TSAdaptCandidatesClear(adapt);CHKERRQ(ierr); 1061 ierr = TSAdaptCandidateAdd(adapt,tab->name,tab->order,1,tab->ccfl,1.*tab->s,PETSC_TRUE);CHKERRQ(ierr); 1062 ierr = TSAdaptChoose(adapt,ts,ts->time_step,&next_scheme,&next_time_step,&accept);CHKERRQ(ierr); 1063 if (accept) { 1064 /* ignore next_scheme for now */ 1065 ts->ptime += ts->time_step; 1066 ts->time_step = next_time_step; 1067 ts->steps++; 1068 ros->status = TS_STEP_COMPLETE; 1069 break; 1070 } else { /* Roll back the current step */ 1071 for (i=0; i<s; i++) w[i] = -tab->bt[i]; 1072 ierr = VecMAXPY(ts->vec_sol,s,w,Y);CHKERRQ(ierr); 1073 ts->time_step = next_time_step; 1074 ros->status = TS_STEP_INCOMPLETE; 1075 } 1076 reject_step: continue; 1077 } 1078 if (ros->status != TS_STEP_COMPLETE && !ts->reason) ts->reason = TS_DIVERGED_STEP_REJECTED; 1079 PetscFunctionReturn(0); 1080 } 1081 1082 #undef __FUNCT__ 1083 #define __FUNCT__ "TSInterpolate_RosW" 1084 static PetscErrorCode TSInterpolate_RosW(TS ts,PetscReal itime,Vec U) 1085 { 1086 TS_RosW *ros = (TS_RosW*)ts->data; 1087 PetscInt s = ros->tableau->s,pinterp = ros->tableau->pinterp,i,j; 1088 PetscReal h; 1089 PetscReal tt,t; 1090 PetscScalar *bt; 1091 const PetscReal *Bt = ros->tableau->binterpt; 1092 PetscErrorCode ierr; 1093 const PetscReal *GammaInv = ros->tableau->GammaInv; 1094 PetscScalar *w = ros->work; 1095 Vec *Y = ros->Y; 1096 1097 PetscFunctionBegin; 1098 if (!Bt) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"TSRosW %s does not have an interpolation formula",ros->tableau->name); 1099 1100 switch (ros->status) { 1101 case TS_STEP_INCOMPLETE: 1102 case TS_STEP_PENDING: 1103 h = ts->time_step; 1104 t = (itime - ts->ptime)/h; 1105 break; 1106 case TS_STEP_COMPLETE: 1107 h = ts->time_step_prev; 1108 t = (itime - ts->ptime)/h + 1; /* In the interval [0,1] */ 1109 break; 1110 default: SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_PLIB,"Invalid TSStepStatus"); 1111 } 1112 ierr = PetscMalloc1(s,&bt);CHKERRQ(ierr); 1113 for (i=0; i<s; i++) bt[i] = 0; 1114 for (j=0,tt=t; j<pinterp; j++,tt*=t) { 1115 for (i=0; i<s; i++) { 1116 bt[i] += Bt[i*pinterp+j] * tt; 1117 } 1118 } 1119 1120 /* y(t+tt*h) = y(t) + Sum bt(tt) * GammaInv * Ydot */ 1121 /*U<-0*/ 1122 ierr = VecZeroEntries(U);CHKERRQ(ierr); 1123 1124 /*U<- Sum bt_i * GammaInv(i,1:i) * Y(1:i) */ 1125 for (j=0; j<s; j++) w[j]=0; 1126 for (j=0; j<s; j++) { 1127 for (i=j; i<s; i++) { 1128 w[j] += bt[i]*GammaInv[i*s+j]; 1129 } 1130 } 1131 ierr = VecMAXPY(U,i,w,Y);CHKERRQ(ierr); 1132 1133 /*X<-y(t) + X*/ 1134 ierr = VecAXPY(U,1.0,ros->VecSolPrev);CHKERRQ(ierr); 1135 1136 ierr = PetscFree(bt);CHKERRQ(ierr); 1137 PetscFunctionReturn(0); 1138 } 1139 1140 /*------------------------------------------------------------*/ 1141 #undef __FUNCT__ 1142 #define __FUNCT__ "TSReset_RosW" 1143 static PetscErrorCode TSReset_RosW(TS ts) 1144 { 1145 TS_RosW *ros = (TS_RosW*)ts->data; 1146 PetscInt s; 1147 PetscErrorCode ierr; 1148 1149 PetscFunctionBegin; 1150 if (!ros->tableau) PetscFunctionReturn(0); 1151 s = ros->tableau->s; 1152 ierr = VecDestroyVecs(s,&ros->Y);CHKERRQ(ierr); 1153 ierr = VecDestroy(&ros->Ydot);CHKERRQ(ierr); 1154 ierr = VecDestroy(&ros->Ystage);CHKERRQ(ierr); 1155 ierr = VecDestroy(&ros->Zdot);CHKERRQ(ierr); 1156 ierr = VecDestroy(&ros->Zstage);CHKERRQ(ierr); 1157 ierr = VecDestroy(&ros->VecSolPrev);CHKERRQ(ierr); 1158 ierr = PetscFree(ros->work);CHKERRQ(ierr); 1159 PetscFunctionReturn(0); 1160 } 1161 1162 #undef __FUNCT__ 1163 #define __FUNCT__ "TSDestroy_RosW" 1164 static PetscErrorCode TSDestroy_RosW(TS ts) 1165 { 1166 PetscErrorCode ierr; 1167 1168 PetscFunctionBegin; 1169 ierr = TSReset_RosW(ts);CHKERRQ(ierr); 1170 ierr = PetscFree(ts->data);CHKERRQ(ierr); 1171 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSRosWGetType_C",NULL);CHKERRQ(ierr); 1172 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSRosWSetType_C",NULL);CHKERRQ(ierr); 1173 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSRosWSetRecomputeJacobian_C",NULL);CHKERRQ(ierr); 1174 PetscFunctionReturn(0); 1175 } 1176 1177 1178 #undef __FUNCT__ 1179 #define __FUNCT__ "TSRosWGetVecs" 1180 static PetscErrorCode TSRosWGetVecs(TS ts,DM dm,Vec *Ydot,Vec *Zdot,Vec *Ystage,Vec *Zstage) 1181 { 1182 TS_RosW *rw = (TS_RosW*)ts->data; 1183 PetscErrorCode ierr; 1184 1185 PetscFunctionBegin; 1186 if (Ydot) { 1187 if (dm && dm != ts->dm) { 1188 ierr = DMGetNamedGlobalVector(dm,"TSRosW_Ydot",Ydot);CHKERRQ(ierr); 1189 } else *Ydot = rw->Ydot; 1190 } 1191 if (Zdot) { 1192 if (dm && dm != ts->dm) { 1193 ierr = DMGetNamedGlobalVector(dm,"TSRosW_Zdot",Zdot);CHKERRQ(ierr); 1194 } else *Zdot = rw->Zdot; 1195 } 1196 if (Ystage) { 1197 if (dm && dm != ts->dm) { 1198 ierr = DMGetNamedGlobalVector(dm,"TSRosW_Ystage",Ystage);CHKERRQ(ierr); 1199 } else *Ystage = rw->Ystage; 1200 } 1201 if (Zstage) { 1202 if (dm && dm != ts->dm) { 1203 ierr = DMGetNamedGlobalVector(dm,"TSRosW_Zstage",Zstage);CHKERRQ(ierr); 1204 } else *Zstage = rw->Zstage; 1205 } 1206 PetscFunctionReturn(0); 1207 } 1208 1209 1210 #undef __FUNCT__ 1211 #define __FUNCT__ "TSRosWRestoreVecs" 1212 static PetscErrorCode TSRosWRestoreVecs(TS ts,DM dm,Vec *Ydot,Vec *Zdot, Vec *Ystage, Vec *Zstage) 1213 { 1214 PetscErrorCode ierr; 1215 1216 PetscFunctionBegin; 1217 if (Ydot) { 1218 if (dm && dm != ts->dm) { 1219 ierr = DMRestoreNamedGlobalVector(dm,"TSRosW_Ydot",Ydot);CHKERRQ(ierr); 1220 } 1221 } 1222 if (Zdot) { 1223 if (dm && dm != ts->dm) { 1224 ierr = DMRestoreNamedGlobalVector(dm,"TSRosW_Zdot",Zdot);CHKERRQ(ierr); 1225 } 1226 } 1227 if (Ystage) { 1228 if (dm && dm != ts->dm) { 1229 ierr = DMRestoreNamedGlobalVector(dm,"TSRosW_Ystage",Ystage);CHKERRQ(ierr); 1230 } 1231 } 1232 if (Zstage) { 1233 if (dm && dm != ts->dm) { 1234 ierr = DMRestoreNamedGlobalVector(dm,"TSRosW_Zstage",Zstage);CHKERRQ(ierr); 1235 } 1236 } 1237 PetscFunctionReturn(0); 1238 } 1239 1240 #undef __FUNCT__ 1241 #define __FUNCT__ "DMCoarsenHook_TSRosW" 1242 static PetscErrorCode DMCoarsenHook_TSRosW(DM fine,DM coarse,void *ctx) 1243 { 1244 PetscFunctionBegin; 1245 PetscFunctionReturn(0); 1246 } 1247 1248 #undef __FUNCT__ 1249 #define __FUNCT__ "DMRestrictHook_TSRosW" 1250 static PetscErrorCode DMRestrictHook_TSRosW(DM fine,Mat restrct,Vec rscale,Mat inject,DM coarse,void *ctx) 1251 { 1252 TS ts = (TS)ctx; 1253 PetscErrorCode ierr; 1254 Vec Ydot,Zdot,Ystage,Zstage; 1255 Vec Ydotc,Zdotc,Ystagec,Zstagec; 1256 1257 PetscFunctionBegin; 1258 ierr = TSRosWGetVecs(ts,fine,&Ydot,&Ystage,&Zdot,&Zstage);CHKERRQ(ierr); 1259 ierr = TSRosWGetVecs(ts,coarse,&Ydotc,&Ystagec,&Zdotc,&Zstagec);CHKERRQ(ierr); 1260 ierr = MatRestrict(restrct,Ydot,Ydotc);CHKERRQ(ierr); 1261 ierr = VecPointwiseMult(Ydotc,rscale,Ydotc);CHKERRQ(ierr); 1262 ierr = MatRestrict(restrct,Ystage,Ystagec);CHKERRQ(ierr); 1263 ierr = VecPointwiseMult(Ystagec,rscale,Ystagec);CHKERRQ(ierr); 1264 ierr = MatRestrict(restrct,Zdot,Zdotc);CHKERRQ(ierr); 1265 ierr = VecPointwiseMult(Zdotc,rscale,Zdotc);CHKERRQ(ierr); 1266 ierr = MatRestrict(restrct,Zstage,Zstagec);CHKERRQ(ierr); 1267 ierr = VecPointwiseMult(Zstagec,rscale,Zstagec);CHKERRQ(ierr); 1268 ierr = TSRosWRestoreVecs(ts,fine,&Ydot,&Ystage,&Zdot,&Zstage);CHKERRQ(ierr); 1269 ierr = TSRosWRestoreVecs(ts,coarse,&Ydotc,&Ystagec,&Zdotc,&Zstagec);CHKERRQ(ierr); 1270 PetscFunctionReturn(0); 1271 } 1272 1273 1274 #undef __FUNCT__ 1275 #define __FUNCT__ "DMSubDomainHook_TSRosW" 1276 static PetscErrorCode DMSubDomainHook_TSRosW(DM fine,DM coarse,void *ctx) 1277 { 1278 PetscFunctionBegin; 1279 PetscFunctionReturn(0); 1280 } 1281 1282 #undef __FUNCT__ 1283 #define __FUNCT__ "DMSubDomainRestrictHook_TSRosW" 1284 static PetscErrorCode DMSubDomainRestrictHook_TSRosW(DM dm,VecScatter gscat,VecScatter lscat,DM subdm,void *ctx) 1285 { 1286 TS ts = (TS)ctx; 1287 PetscErrorCode ierr; 1288 Vec Ydot,Zdot,Ystage,Zstage; 1289 Vec Ydots,Zdots,Ystages,Zstages; 1290 1291 PetscFunctionBegin; 1292 ierr = TSRosWGetVecs(ts,dm,&Ydot,&Ystage,&Zdot,&Zstage);CHKERRQ(ierr); 1293 ierr = TSRosWGetVecs(ts,subdm,&Ydots,&Ystages,&Zdots,&Zstages);CHKERRQ(ierr); 1294 1295 ierr = VecScatterBegin(gscat,Ydot,Ydots,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1296 ierr = VecScatterEnd(gscat,Ydot,Ydots,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1297 1298 ierr = VecScatterBegin(gscat,Ystage,Ystages,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1299 ierr = VecScatterEnd(gscat,Ystage,Ystages,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1300 1301 ierr = VecScatterBegin(gscat,Zdot,Zdots,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1302 ierr = VecScatterEnd(gscat,Zdot,Zdots,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1303 1304 ierr = VecScatterBegin(gscat,Zstage,Zstages,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1305 ierr = VecScatterEnd(gscat,Zstage,Zstages,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1306 1307 ierr = TSRosWRestoreVecs(ts,dm,&Ydot,&Ystage,&Zdot,&Zstage);CHKERRQ(ierr); 1308 ierr = TSRosWRestoreVecs(ts,subdm,&Ydots,&Ystages,&Zdots,&Zstages);CHKERRQ(ierr); 1309 PetscFunctionReturn(0); 1310 } 1311 1312 /* 1313 This defines the nonlinear equation that is to be solved with SNES 1314 G(U) = F[t0+Theta*dt, U, (U-U0)*shift] = 0 1315 */ 1316 #undef __FUNCT__ 1317 #define __FUNCT__ "SNESTSFormFunction_RosW" 1318 static PetscErrorCode SNESTSFormFunction_RosW(SNES snes,Vec U,Vec F,TS ts) 1319 { 1320 TS_RosW *ros = (TS_RosW*)ts->data; 1321 PetscErrorCode ierr; 1322 Vec Ydot,Zdot,Ystage,Zstage; 1323 PetscReal shift = ros->scoeff / ts->time_step; 1324 DM dm,dmsave; 1325 1326 PetscFunctionBegin; 1327 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 1328 ierr = TSRosWGetVecs(ts,dm,&Ydot,&Zdot,&Ystage,&Zstage);CHKERRQ(ierr); 1329 ierr = VecWAXPY(Ydot,shift,U,Zdot);CHKERRQ(ierr); /* Ydot = shift*U + Zdot */ 1330 ierr = VecWAXPY(Ystage,1.0,U,Zstage);CHKERRQ(ierr); /* Ystage = U + Zstage */ 1331 dmsave = ts->dm; 1332 ts->dm = dm; 1333 ierr = TSComputeIFunction(ts,ros->stage_time,Ystage,Ydot,F,PETSC_FALSE);CHKERRQ(ierr); 1334 ts->dm = dmsave; 1335 ierr = TSRosWRestoreVecs(ts,dm,&Ydot,&Zdot,&Ystage,&Zstage);CHKERRQ(ierr); 1336 PetscFunctionReturn(0); 1337 } 1338 1339 #undef __FUNCT__ 1340 #define __FUNCT__ "SNESTSFormJacobian_RosW" 1341 static PetscErrorCode SNESTSFormJacobian_RosW(SNES snes,Vec U,Mat *A,Mat *B,MatStructure *str,TS ts) 1342 { 1343 TS_RosW *ros = (TS_RosW*)ts->data; 1344 Vec Ydot,Zdot,Ystage,Zstage; 1345 PetscReal shift = ros->scoeff / ts->time_step; 1346 PetscErrorCode ierr; 1347 DM dm,dmsave; 1348 1349 PetscFunctionBegin; 1350 /* ros->Ydot and ros->Ystage have already been computed in SNESTSFormFunction_RosW (SNES guarantees this) */ 1351 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 1352 ierr = TSRosWGetVecs(ts,dm,&Ydot,&Zdot,&Ystage,&Zstage);CHKERRQ(ierr); 1353 dmsave = ts->dm; 1354 ts->dm = dm; 1355 ierr = TSComputeIJacobian(ts,ros->stage_time,Ystage,Ydot,shift,A,B,str,PETSC_TRUE);CHKERRQ(ierr); 1356 ts->dm = dmsave; 1357 ierr = TSRosWRestoreVecs(ts,dm,&Ydot,&Zdot,&Ystage,&Zstage);CHKERRQ(ierr); 1358 PetscFunctionReturn(0); 1359 } 1360 1361 #undef __FUNCT__ 1362 #define __FUNCT__ "TSSetUp_RosW" 1363 static PetscErrorCode TSSetUp_RosW(TS ts) 1364 { 1365 TS_RosW *ros = (TS_RosW*)ts->data; 1366 RosWTableau tab = ros->tableau; 1367 PetscInt s = tab->s; 1368 PetscErrorCode ierr; 1369 DM dm; 1370 1371 PetscFunctionBegin; 1372 if (!ros->tableau) { 1373 ierr = TSRosWSetType(ts,TSRosWDefault);CHKERRQ(ierr); 1374 } 1375 ierr = VecDuplicateVecs(ts->vec_sol,s,&ros->Y);CHKERRQ(ierr); 1376 ierr = VecDuplicate(ts->vec_sol,&ros->Ydot);CHKERRQ(ierr); 1377 ierr = VecDuplicate(ts->vec_sol,&ros->Ystage);CHKERRQ(ierr); 1378 ierr = VecDuplicate(ts->vec_sol,&ros->Zdot);CHKERRQ(ierr); 1379 ierr = VecDuplicate(ts->vec_sol,&ros->Zstage);CHKERRQ(ierr); 1380 ierr = VecDuplicate(ts->vec_sol,&ros->VecSolPrev);CHKERRQ(ierr); 1381 ierr = PetscMalloc1(s,&ros->work);CHKERRQ(ierr); 1382 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 1383 if (dm) { 1384 ierr = DMCoarsenHookAdd(dm,DMCoarsenHook_TSRosW,DMRestrictHook_TSRosW,ts);CHKERRQ(ierr); 1385 ierr = DMSubDomainHookAdd(dm,DMSubDomainHook_TSRosW,DMSubDomainRestrictHook_TSRosW,ts);CHKERRQ(ierr); 1386 } 1387 PetscFunctionReturn(0); 1388 } 1389 /*------------------------------------------------------------*/ 1390 1391 #undef __FUNCT__ 1392 #define __FUNCT__ "TSSetFromOptions_RosW" 1393 static PetscErrorCode TSSetFromOptions_RosW(TS ts) 1394 { 1395 TS_RosW *ros = (TS_RosW*)ts->data; 1396 PetscErrorCode ierr; 1397 char rostype[256]; 1398 1399 PetscFunctionBegin; 1400 ierr = PetscOptionsHead("RosW ODE solver options");CHKERRQ(ierr); 1401 { 1402 RosWTableauLink link; 1403 PetscInt count,choice; 1404 PetscBool flg; 1405 const char **namelist; 1406 SNES snes; 1407 1408 ierr = PetscStrncpy(rostype,TSRosWDefault,sizeof(rostype));CHKERRQ(ierr); 1409 for (link=RosWTableauList,count=0; link; link=link->next,count++) ; 1410 ierr = PetscMalloc1(count,&namelist);CHKERRQ(ierr); 1411 for (link=RosWTableauList,count=0; link; link=link->next,count++) namelist[count] = link->tab.name; 1412 ierr = PetscOptionsEList("-ts_rosw_type","Family of Rosenbrock-W method","TSRosWSetType",(const char*const*)namelist,count,rostype,&choice,&flg);CHKERRQ(ierr); 1413 ierr = TSRosWSetType(ts,flg ? namelist[choice] : rostype);CHKERRQ(ierr); 1414 ierr = PetscFree(namelist);CHKERRQ(ierr); 1415 1416 ierr = PetscOptionsBool("-ts_rosw_recompute_jacobian","Recompute the Jacobian at each stage","TSRosWSetRecomputeJacobian",ros->recompute_jacobian,&ros->recompute_jacobian,NULL);CHKERRQ(ierr); 1417 1418 /* Rosenbrock methods are linearly implicit, so set that unless the user has specifically asked for something else */ 1419 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 1420 if (!((PetscObject)snes)->type_name) { 1421 ierr = SNESSetType(snes,SNESKSPONLY);CHKERRQ(ierr); 1422 } 1423 ierr = SNESSetFromOptions(snes);CHKERRQ(ierr); 1424 } 1425 ierr = PetscOptionsTail();CHKERRQ(ierr); 1426 PetscFunctionReturn(0); 1427 } 1428 1429 #undef __FUNCT__ 1430 #define __FUNCT__ "PetscFormatRealArray" 1431 static PetscErrorCode PetscFormatRealArray(char buf[],size_t len,const char *fmt,PetscInt n,const PetscReal x[]) 1432 { 1433 PetscErrorCode ierr; 1434 PetscInt i; 1435 size_t left,count; 1436 char *p; 1437 1438 PetscFunctionBegin; 1439 for (i=0,p=buf,left=len; i<n; i++) { 1440 ierr = PetscSNPrintfCount(p,left,fmt,&count,x[i]);CHKERRQ(ierr); 1441 if (count >= left) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Insufficient space in buffer"); 1442 left -= count; 1443 p += count; 1444 *p++ = ' '; 1445 } 1446 p[i ? 0 : -1] = 0; 1447 PetscFunctionReturn(0); 1448 } 1449 1450 #undef __FUNCT__ 1451 #define __FUNCT__ "TSView_RosW" 1452 static PetscErrorCode TSView_RosW(TS ts,PetscViewer viewer) 1453 { 1454 TS_RosW *ros = (TS_RosW*)ts->data; 1455 RosWTableau tab = ros->tableau; 1456 PetscBool iascii; 1457 PetscErrorCode ierr; 1458 TSAdapt adapt; 1459 1460 PetscFunctionBegin; 1461 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 1462 if (iascii) { 1463 TSRosWType rostype; 1464 PetscInt i; 1465 PetscReal abscissa[512]; 1466 char buf[512]; 1467 ierr = TSRosWGetType(ts,&rostype);CHKERRQ(ierr); 1468 ierr = PetscViewerASCIIPrintf(viewer," Rosenbrock-W %s\n",rostype);CHKERRQ(ierr); 1469 ierr = PetscFormatRealArray(buf,sizeof(buf),"% 8.6f",tab->s,tab->ASum);CHKERRQ(ierr); 1470 ierr = PetscViewerASCIIPrintf(viewer," Abscissa of A = %s\n",buf);CHKERRQ(ierr); 1471 for (i=0; i<tab->s; i++) abscissa[i] = tab->ASum[i] + tab->Gamma[i]; 1472 ierr = PetscFormatRealArray(buf,sizeof(buf),"% 8.6f",tab->s,abscissa);CHKERRQ(ierr); 1473 ierr = PetscViewerASCIIPrintf(viewer," Abscissa of A+Gamma = %s\n",buf);CHKERRQ(ierr); 1474 } 1475 ierr = TSGetAdapt(ts,&adapt);CHKERRQ(ierr); 1476 ierr = TSAdaptView(adapt,viewer);CHKERRQ(ierr); 1477 ierr = SNESView(ts->snes,viewer);CHKERRQ(ierr); 1478 PetscFunctionReturn(0); 1479 } 1480 1481 #undef __FUNCT__ 1482 #define __FUNCT__ "TSLoad_RosW" 1483 static PetscErrorCode TSLoad_RosW(TS ts,PetscViewer viewer) 1484 { 1485 PetscErrorCode ierr; 1486 SNES snes; 1487 TSAdapt tsadapt; 1488 1489 PetscFunctionBegin; 1490 ierr = TSGetAdapt(ts,&tsadapt);CHKERRQ(ierr); 1491 ierr = TSAdaptLoad(tsadapt,viewer);CHKERRQ(ierr); 1492 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 1493 ierr = SNESLoad(snes,viewer);CHKERRQ(ierr); 1494 /* function and Jacobian context for SNES when used with TS is always ts object */ 1495 ierr = SNESSetFunction(snes,NULL,NULL,ts);CHKERRQ(ierr); 1496 ierr = SNESSetJacobian(snes,NULL,NULL,NULL,ts);CHKERRQ(ierr); 1497 PetscFunctionReturn(0); 1498 } 1499 1500 #undef __FUNCT__ 1501 #define __FUNCT__ "TSRosWSetType" 1502 /*@C 1503 TSRosWSetType - Set the type of Rosenbrock-W scheme 1504 1505 Logically collective 1506 1507 Input Parameter: 1508 + ts - timestepping context 1509 - rostype - type of Rosenbrock-W scheme 1510 1511 Level: beginner 1512 1513 .seealso: TSRosWGetType(), TSROSW, TSROSW2M, TSROSW2P, TSROSWRA3PW, TSROSWRA34PW2, TSROSWRODAS3, TSROSWSANDU3, TSROSWASSP3P3S1C, TSROSWLASSP3P4S2C, TSROSWLLSSP3P4S2C, TSROSWARK3 1514 @*/ 1515 PetscErrorCode TSRosWSetType(TS ts,TSRosWType rostype) 1516 { 1517 PetscErrorCode ierr; 1518 1519 PetscFunctionBegin; 1520 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1521 ierr = PetscTryMethod(ts,"TSRosWSetType_C",(TS,TSRosWType),(ts,rostype));CHKERRQ(ierr); 1522 PetscFunctionReturn(0); 1523 } 1524 1525 #undef __FUNCT__ 1526 #define __FUNCT__ "TSRosWGetType" 1527 /*@C 1528 TSRosWGetType - Get the type of Rosenbrock-W scheme 1529 1530 Logically collective 1531 1532 Input Parameter: 1533 . ts - timestepping context 1534 1535 Output Parameter: 1536 . rostype - type of Rosenbrock-W scheme 1537 1538 Level: intermediate 1539 1540 .seealso: TSRosWGetType() 1541 @*/ 1542 PetscErrorCode TSRosWGetType(TS ts,TSRosWType *rostype) 1543 { 1544 PetscErrorCode ierr; 1545 1546 PetscFunctionBegin; 1547 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1548 ierr = PetscUseMethod(ts,"TSRosWGetType_C",(TS,TSRosWType*),(ts,rostype));CHKERRQ(ierr); 1549 PetscFunctionReturn(0); 1550 } 1551 1552 #undef __FUNCT__ 1553 #define __FUNCT__ "TSRosWSetRecomputeJacobian" 1554 /*@C 1555 TSRosWSetRecomputeJacobian - Set whether to recompute the Jacobian at each stage. The default is to update the Jacobian once per step. 1556 1557 Logically collective 1558 1559 Input Parameter: 1560 + ts - timestepping context 1561 - flg - PETSC_TRUE to recompute the Jacobian at each stage 1562 1563 Level: intermediate 1564 1565 .seealso: TSRosWGetType() 1566 @*/ 1567 PetscErrorCode TSRosWSetRecomputeJacobian(TS ts,PetscBool flg) 1568 { 1569 PetscErrorCode ierr; 1570 1571 PetscFunctionBegin; 1572 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1573 ierr = PetscTryMethod(ts,"TSRosWSetRecomputeJacobian_C",(TS,PetscBool),(ts,flg));CHKERRQ(ierr); 1574 PetscFunctionReturn(0); 1575 } 1576 1577 #undef __FUNCT__ 1578 #define __FUNCT__ "TSRosWGetType_RosW" 1579 PetscErrorCode TSRosWGetType_RosW(TS ts,TSRosWType *rostype) 1580 { 1581 TS_RosW *ros = (TS_RosW*)ts->data; 1582 PetscErrorCode ierr; 1583 1584 PetscFunctionBegin; 1585 if (!ros->tableau) {ierr = TSRosWSetType(ts,TSRosWDefault);CHKERRQ(ierr);} 1586 *rostype = ros->tableau->name; 1587 PetscFunctionReturn(0); 1588 } 1589 1590 #undef __FUNCT__ 1591 #define __FUNCT__ "TSRosWSetType_RosW" 1592 PetscErrorCode TSRosWSetType_RosW(TS ts,TSRosWType rostype) 1593 { 1594 TS_RosW *ros = (TS_RosW*)ts->data; 1595 PetscErrorCode ierr; 1596 PetscBool match; 1597 RosWTableauLink link; 1598 1599 PetscFunctionBegin; 1600 if (ros->tableau) { 1601 ierr = PetscStrcmp(ros->tableau->name,rostype,&match);CHKERRQ(ierr); 1602 if (match) PetscFunctionReturn(0); 1603 } 1604 for (link = RosWTableauList; link; link=link->next) { 1605 ierr = PetscStrcmp(link->tab.name,rostype,&match);CHKERRQ(ierr); 1606 if (match) { 1607 ierr = TSReset_RosW(ts);CHKERRQ(ierr); 1608 ros->tableau = &link->tab; 1609 PetscFunctionReturn(0); 1610 } 1611 } 1612 SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_UNKNOWN_TYPE,"Could not find '%s'",rostype); 1613 PetscFunctionReturn(0); 1614 } 1615 1616 #undef __FUNCT__ 1617 #define __FUNCT__ "TSRosWSetRecomputeJacobian_RosW" 1618 PetscErrorCode TSRosWSetRecomputeJacobian_RosW(TS ts,PetscBool flg) 1619 { 1620 TS_RosW *ros = (TS_RosW*)ts->data; 1621 1622 PetscFunctionBegin; 1623 ros->recompute_jacobian = flg; 1624 PetscFunctionReturn(0); 1625 } 1626 1627 1628 /* ------------------------------------------------------------ */ 1629 /*MC 1630 TSROSW - ODE solver using Rosenbrock-W schemes 1631 1632 These methods are intended for problems with well-separated time scales, especially when a slow scale is strongly 1633 nonlinear such that it is expensive to solve with a fully implicit method. The user should provide the stiff part 1634 of the equation using TSSetIFunction() and the non-stiff part with TSSetRHSFunction(). 1635 1636 Notes: 1637 This method currently only works with autonomous ODE and DAE. 1638 1639 Developer notes: 1640 Rosenbrock-W methods are typically specified for autonomous ODE 1641 1642 $ udot = f(u) 1643 1644 by the stage equations 1645 1646 $ k_i = h f(u_0 + sum_j alpha_ij k_j) + h J sum_j gamma_ij k_j 1647 1648 and step completion formula 1649 1650 $ u_1 = u_0 + sum_j b_j k_j 1651 1652 with step size h and coefficients alpha_ij, gamma_ij, and b_i. Implementing the method in this form would require f(u) 1653 and the Jacobian J to be available, in addition to the shifted matrix I - h gamma_ii J. Following Hairer and Wanner, 1654 we define new variables for the stage equations 1655 1656 $ y_i = gamma_ij k_j 1657 1658 The k_j can be recovered because Gamma is invertible. Let C be the lower triangular part of Gamma^{-1} and define 1659 1660 $ A = Alpha Gamma^{-1}, bt^T = b^T Gamma^{-i} 1661 1662 to rewrite the method as 1663 1664 $ [M/(h gamma_ii) - J] y_i = f(u_0 + sum_j a_ij y_j) + M sum_j (c_ij/h) y_j 1665 $ u_1 = u_0 + sum_j bt_j y_j 1666 1667 where we have introduced the mass matrix M. Continue by defining 1668 1669 $ ydot_i = 1/(h gamma_ii) y_i - sum_j (c_ij/h) y_j 1670 1671 or, more compactly in tensor notation 1672 1673 $ Ydot = 1/h (Gamma^{-1} \otimes I) Y . 1674 1675 Note that Gamma^{-1} is lower triangular. With this definition of Ydot in terms of known quantities and the current 1676 stage y_i, the stage equations reduce to performing one Newton step (typically with a lagged Jacobian) on the 1677 equation 1678 1679 $ g(u_0 + sum_j a_ij y_j + y_i, ydot_i) = 0 1680 1681 with initial guess y_i = 0. 1682 1683 Level: beginner 1684 1685 .seealso: TSCreate(), TS, TSSetType(), TSRosWSetType(), TSRosWRegister(), TSROSW2M, TSROSW2P, TSROSWRA3PW, TSROSWRA34PW2, TSROSWRODAS3, 1686 TSROSWSANDU3, TSROSWASSP3P3S1C, TSROSWLASSP3P4S2C, TSROSWLLSSP3P4S2C, TSROSWGRK4T, TSROSWSHAMP4, TSROSWVELDD4, TSROSW4L 1687 M*/ 1688 #undef __FUNCT__ 1689 #define __FUNCT__ "TSCreate_RosW" 1690 PETSC_EXTERN PetscErrorCode TSCreate_RosW(TS ts) 1691 { 1692 TS_RosW *ros; 1693 PetscErrorCode ierr; 1694 1695 PetscFunctionBegin; 1696 ierr = TSRosWInitializePackage();CHKERRQ(ierr); 1697 1698 ts->ops->reset = TSReset_RosW; 1699 ts->ops->destroy = TSDestroy_RosW; 1700 ts->ops->view = TSView_RosW; 1701 ts->ops->load = TSLoad_RosW; 1702 ts->ops->setup = TSSetUp_RosW; 1703 ts->ops->step = TSStep_RosW; 1704 ts->ops->interpolate = TSInterpolate_RosW; 1705 ts->ops->evaluatestep = TSEvaluateStep_RosW; 1706 ts->ops->setfromoptions = TSSetFromOptions_RosW; 1707 ts->ops->snesfunction = SNESTSFormFunction_RosW; 1708 ts->ops->snesjacobian = SNESTSFormJacobian_RosW; 1709 1710 ierr = PetscNewLog(ts,TS_RosW,&ros);CHKERRQ(ierr); 1711 ts->data = (void*)ros; 1712 1713 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSRosWGetType_C",TSRosWGetType_RosW);CHKERRQ(ierr); 1714 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSRosWSetType_C",TSRosWSetType_RosW);CHKERRQ(ierr); 1715 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSRosWSetRecomputeJacobian_C",TSRosWSetRecomputeJacobian_RosW);CHKERRQ(ierr); 1716 PetscFunctionReturn(0); 1717 } 1718