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