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