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