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