1 /* 2 Code for timestepping with Rosenbrock W methods 3 4 Notes: 5 The general system is written as 6 7 G(t,X,Xdot) = F(t,X) 8 9 where G represents the stiff part of the physics and F represents the non-stiff part. 10 This method is designed to be linearly implicit on G and can use an approximate and lagged Jacobian. 11 12 */ 13 #include <petsc-private/tsimpl.h> /*I "petscts.h" I*/ 14 15 #include <../src/mat/blockinvert.h> 16 17 static const TSRosWType TSRosWDefault = TSROSWRA34PW2; 18 static PetscBool TSRosWRegisterAllCalled; 19 static PetscBool TSRosWPackageInitialized; 20 21 typedef struct _RosWTableau *RosWTableau; 22 struct _RosWTableau { 23 char *name; 24 PetscInt order; /* Classical approximation order of the method */ 25 PetscInt s; /* Number of stages */ 26 PetscInt pinterp; /* Interpolation order */ 27 PetscReal *A; /* Propagation table, strictly lower triangular */ 28 PetscReal *Gamma; /* Stage table, lower triangular with nonzero diagonal */ 29 PetscBool *GammaZeroDiag; /* Diagonal entries that are zero in stage table Gamma, vector indicating explicit statages */ 30 PetscReal *GammaExplicitCorr; /* Coefficients for correction terms needed for explicit stages in transformed variables*/ 31 PetscReal *b; /* Step completion table */ 32 PetscReal *bembed; /* Step completion table for embedded method of order one less */ 33 PetscReal *ASum; /* Row sum of A */ 34 PetscReal *GammaSum; /* Row sum of Gamma, only needed for non-autonomous systems */ 35 PetscReal *At; /* Propagation table in transformed variables */ 36 PetscReal *bt; /* Step completion table in transformed variables */ 37 PetscReal *bembedt; /* Step completion table of order one less in transformed variables */ 38 PetscReal *GammaInv; /* Inverse of Gamma, used for transformed variables */ 39 PetscReal ccfl; /* Placeholder for CFL coefficient relative to forward Euler */ 40 PetscReal *binterpt; /* Dense output formula */ 41 }; 42 typedef struct _RosWTableauLink *RosWTableauLink; 43 struct _RosWTableauLink { 44 struct _RosWTableau tab; 45 RosWTableauLink next; 46 }; 47 static RosWTableauLink RosWTableauList; 48 49 typedef struct { 50 RosWTableau tableau; 51 Vec *Y; /* States computed during the step, used to complete the step */ 52 Vec Ydot; /* Work vector holding Ydot during residual evaluation */ 53 Vec Ystage; /* Work vector for the state value at each stage */ 54 Vec Zdot; /* Ydot = Zdot + shift*Y */ 55 Vec Zstage; /* Y = Zstage + Y */ 56 Vec VecSolPrev; /* Work vector holding the solution from the previous step (used for interpolation)*/ 57 PetscScalar *work; /* Scalar work space of length number of stages, used to prepare VecMAXPY() */ 58 PetscReal shift; 59 PetscReal stage_time; 60 PetscReal stage_explicit; /* Flag indicates that the current stage is explicit */ 61 PetscBool recompute_jacobian; /* Recompute the Jacobian at each stage, default is to freeze the Jacobian at the start of each step */ 62 TSStepStatus status; 63 } TS_RosW; 64 65 /*MC 66 TSROSWTHETA1 - One stage first order L-stable Rosenbrock-W scheme (aka theta method). 67 68 Only an approximate Jacobian is needed. 69 70 Level: intermediate 71 72 .seealso: TSROSW 73 M*/ 74 75 /*MC 76 TSROSWTHETA2 - One stage second order A-stable Rosenbrock-W scheme (aka theta method). 77 78 Only an approximate Jacobian is needed. 79 80 Level: intermediate 81 82 .seealso: TSROSW 83 M*/ 84 85 /*MC 86 TSROSW2M - Two stage second order L-stable Rosenbrock-W scheme. 87 88 Only an approximate Jacobian is needed. By default, it is only recomputed once per step. This method is a reflection of TSROSW2P. 89 90 Level: intermediate 91 92 .seealso: TSROSW 93 M*/ 94 95 /*MC 96 TSROSW2P - Two stage second order L-stable Rosenbrock-W scheme. 97 98 Only an approximate Jacobian is needed. By default, it is only recomputed once per step. This method is a reflection of TSROSW2M. 99 100 Level: intermediate 101 102 .seealso: TSROSW 103 M*/ 104 105 /*MC 106 TSROSWRA3PW - Three stage third order Rosenbrock-W scheme for PDAE of index 1. 107 108 Only an approximate Jacobian is needed. By default, it is only recomputed once per step. 109 110 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. 111 112 References: 113 Rang and Angermann, New Rosenbrock-W methods of order 3 for partial differential algebraic equations of index 1, 2005. 114 115 Level: intermediate 116 117 .seealso: TSROSW 118 M*/ 119 120 /*MC 121 TSROSWRA34PW2 - Four stage third order L-stable Rosenbrock-W scheme for PDAE of index 1. 122 123 Only an approximate Jacobian is needed. By default, it is only recomputed once per step. 124 125 This is strongly A-stable with R(infty) = 0. The embedded method of order 2 is strongly A-stable with R(infty) = 0.48. 126 127 References: 128 Rang and Angermann, New Rosenbrock-W methods of order 3 for partial differential algebraic equations of index 1, 2005. 129 130 Level: intermediate 131 132 .seealso: TSROSW 133 M*/ 134 135 /*MC 136 TSROSWRODAS3 - Four stage third order L-stable Rosenbrock scheme 137 138 By default, the Jacobian is only recomputed once per step. 139 140 Both the third order and embedded second order methods are stiffly accurate and L-stable. 141 142 References: 143 Sandu et al, Benchmarking stiff ODE solvers for atmospheric chemistry problems II, Rosenbrock solvers, 1997. 144 145 Level: intermediate 146 147 .seealso: TSROSW, TSROSWSANDU3 148 M*/ 149 150 /*MC 151 TSROSWSANDU3 - Three stage third order L-stable Rosenbrock scheme 152 153 By default, the Jacobian is only recomputed once per step. 154 155 The third order method is L-stable, but not stiffly accurate. 156 The second order embedded method is strongly A-stable with R(infty) = 0.5. 157 The internal stages are L-stable. 158 This method is called ROS3 in the paper. 159 160 References: 161 Sandu et al, Benchmarking stiff ODE solvers for atmospheric chemistry problems II, Rosenbrock solvers, 1997. 162 163 Level: intermediate 164 165 .seealso: TSROSW, TSROSWRODAS3 166 M*/ 167 168 /*MC 169 TSROSWASSP3P3S1C - A-stable Rosenbrock-W method with SSP explicit part, third order, three stages 170 171 By default, the Jacobian is only recomputed once per step. 172 173 A-stable SPP explicit order 3, 3 stages, CFL 1 (eff = 1/3) 174 175 References: 176 Emil Constantinescu 177 178 Level: intermediate 179 180 .seealso: TSROSW, TSROSWLASSP3P4S2C, TSROSWLLSSP3P4S2C, SSP 181 M*/ 182 183 /*MC 184 TSROSWLASSP3P4S2C - L-stable Rosenbrock-W method with SSP explicit part, third order, four stages 185 186 By default, the Jacobian is only recomputed once per step. 187 188 L-stable (A-stable embedded) SPP explicit order 3, 4 stages, CFL 2 (eff = 1/2) 189 190 References: 191 Emil Constantinescu 192 193 Level: intermediate 194 195 .seealso: TSROSW, TSROSWASSP3P3S1C, TSROSWLLSSP3P4S2C, TSSSP 196 M*/ 197 198 /*MC 199 TSROSWLLSSP3P4S2C - L-stable Rosenbrock-W method with SSP explicit part, third order, four stages 200 201 By default, the Jacobian is only recomputed once per step. 202 203 L-stable (L-stable embedded) SPP explicit order 3, 4 stages, CFL 2 (eff = 1/2) 204 205 References: 206 Emil Constantinescu 207 208 Level: intermediate 209 210 .seealso: TSROSW, TSROSWASSP3P3S1C, TSROSWLASSP3P4S2C, TSSSP 211 M*/ 212 213 #undef __FUNCT__ 214 #define __FUNCT__ "TSRosWRegisterAll" 215 /*@C 216 TSRosWRegisterAll - Registers all of the additive Runge-Kutta implicit-explicit methods in TSRosW 217 218 Not Collective, but should be called by all processes which will need the schemes to be registered 219 220 Level: advanced 221 222 .keywords: TS, TSRosW, register, all 223 224 .seealso: TSRosWRegisterDestroy() 225 @*/ 226 PetscErrorCode TSRosWRegisterAll(void) 227 { 228 PetscErrorCode ierr; 229 230 PetscFunctionBegin; 231 if (TSRosWRegisterAllCalled) PetscFunctionReturn(0); 232 TSRosWRegisterAllCalled = PETSC_TRUE; 233 234 { 235 const PetscReal 236 A = 0, 237 Gamma = 1, 238 b = 1, 239 binterpt=1; 240 241 ierr = TSRosWRegister(TSROSWTHETA1,1,1,&A,&Gamma,&b,PETSC_NULL,1,&binterpt);CHKERRQ(ierr); 242 } 243 244 { 245 const PetscReal 246 A= 0, 247 Gamma = 0.5, 248 b = 1, 249 binterpt=1; 250 ierr = TSRosWRegister(TSROSWTHETA2,2,1,&A,&Gamma,&b,PETSC_NULL,1,&binterpt);CHKERRQ(ierr); 251 } 252 253 { 254 const PetscReal g = 1. + 1./PetscSqrtReal(2.0); 255 const PetscReal 256 A[2][2] = {{0,0}, {1.,0}}, 257 Gamma[2][2] = {{g,0}, {-2.*g,g}}, 258 b[2] = {0.5,0.5}, 259 b1[2] = {1.0,0.0}; 260 PetscReal binterpt[2][2]; 261 binterpt[0][0]=g-1.0; 262 binterpt[1][0]=2.0-g; 263 binterpt[0][1]=g-1.5; 264 binterpt[1][1]=1.5-g; 265 ierr = TSRosWRegister(TSROSW2P,2,2,&A[0][0],&Gamma[0][0],b,b1,2,&binterpt[0][0]);CHKERRQ(ierr); 266 } 267 { 268 const PetscReal g = 1. - 1./PetscSqrtReal(2.0); 269 const PetscReal 270 A[2][2] = {{0,0}, {1.,0}}, 271 Gamma[2][2] = {{g,0}, {-2.*g,g}}, 272 b[2] = {0.5,0.5}, 273 b1[2] = {1.0,0.0}; 274 PetscReal binterpt[2][2]; 275 binterpt[0][0]=g-1.0; 276 binterpt[1][0]=2.0-g; 277 binterpt[0][1]=g-1.5; 278 binterpt[1][1]=1.5-g; 279 ierr = TSRosWRegister(TSROSW2M,2,2,&A[0][0],&Gamma[0][0],b,b1,2,&binterpt[0][0]);CHKERRQ(ierr); 280 } 281 { 282 const PetscReal g = 7.8867513459481287e-01; 283 PetscReal binterpt[3][2]; 284 const PetscReal 285 A[3][3] = {{0,0,0}, 286 {1.5773502691896257e+00,0,0}, 287 {0.5,0,0}}, 288 Gamma[3][3] = {{g,0,0}, 289 {-1.5773502691896257e+00,g,0}, 290 {-6.7075317547305480e-01,-1.7075317547305482e-01,g}}, 291 b[3] = {1.0566243270259355e-01,4.9038105676657971e-02,8.4529946162074843e-01}, 292 b2[3] = {-1.7863279495408180e-01,1./3.,8.4529946162074843e-01}; 293 294 binterpt[0][0]=-0.8094010767585034; 295 binterpt[1][0]=-0.5; 296 binterpt[2][0]=2.3094010767585034; 297 binterpt[0][1]=0.9641016151377548; 298 binterpt[1][1]=0.5; 299 binterpt[2][1]=-1.4641016151377548; 300 ierr = TSRosWRegister(TSROSWRA3PW,3,3,&A[0][0],&Gamma[0][0],b,b2,2,&binterpt[0][0]);CHKERRQ(ierr); 301 } 302 { 303 PetscReal binterpt[4][3]; 304 const PetscReal g = 4.3586652150845900e-01; 305 const PetscReal 306 A[4][4] = {{0,0,0,0}, 307 {8.7173304301691801e-01,0,0,0}, 308 {8.4457060015369423e-01,-1.1299064236484185e-01,0,0}, 309 {0,0,1.,0}}, 310 Gamma[4][4] = {{g,0,0,0}, 311 {-8.7173304301691801e-01,g,0,0}, 312 {-9.0338057013044082e-01,5.4180672388095326e-02,g,0}, 313 {2.4212380706095346e-01,-1.2232505839045147e+00,5.4526025533510214e-01,g}}, 314 b[4] = {2.4212380706095346e-01,-1.2232505839045147e+00,1.5452602553351020e+00,4.3586652150845900e-01}, 315 b2[4] = {3.7810903145819369e-01,-9.6042292212423178e-02,5.0000000000000000e-01,2.1793326075422950e-01}; 316 317 binterpt[0][0]=1.0564298455794094; 318 binterpt[1][0]=2.296429974281067; 319 binterpt[2][0]=-1.307599564525376; 320 binterpt[3][0]=-1.045260255335102; 321 binterpt[0][1]=-1.3864882699759573; 322 binterpt[1][1]=-8.262611700275677; 323 binterpt[2][1]=7.250979895056055; 324 binterpt[3][1]=2.398120075195581; 325 binterpt[0][2]=0.5721822314575016; 326 binterpt[1][2]=4.742931142090097; 327 binterpt[2][2]=-4.398120075195578; 328 binterpt[3][2]=-0.9169932983520199; 329 330 ierr = TSRosWRegister(TSROSWRA34PW2,3,4,&A[0][0],&Gamma[0][0],b,b2,3,&binterpt[0][0]);CHKERRQ(ierr); 331 } 332 { 333 const PetscReal g = 0.5; 334 const PetscReal 335 A[4][4] = {{0,0,0,0}, 336 {0,0,0,0}, 337 {1.,0,0,0}, 338 {0.75,-0.25,0.5,0}}, 339 Gamma[4][4] = {{g,0,0,0}, 340 {1.,g,0,0}, 341 {-0.25,-0.25,g,0}, 342 {1./12,1./12,-2./3,g}}, 343 b[4] = {5./6,-1./6,-1./6,0.5}, 344 b2[4] = {0.75,-0.25,0.5,0}; 345 ierr = TSRosWRegister(TSROSWRODAS3,3,4,&A[0][0],&Gamma[0][0],b,b2,0,PETSC_NULL);CHKERRQ(ierr); 346 } 347 { 348 const PetscReal g = 0.43586652150845899941601945119356; 349 const PetscReal 350 A[3][3] = {{0,0,0}, 351 {g,0,0}, 352 {g,0,0}}, 353 Gamma[3][3] = {{g,0,0}, 354 {-0.19294655696029095575009695436041,g,0}, 355 {0,1.74927148125794685173529749738960,g}}, 356 b[3] = {-0.75457412385404315829818998646589,1.94100407061964420292840123379419,-0.18642994676560104463021124732829}, 357 b2[3] = {-1.53358745784149585370766523913002,2.81745131148625772213931745457622,-0.28386385364476186843165221544619}; 358 359 PetscReal binterpt[3][2]; 360 binterpt[0][0]=3.793692883777660870425141387941; 361 binterpt[1][0]=-2.918692883777660870425141387941; 362 binterpt[2][0]=0.125; 363 binterpt[0][1]=-0.725741064379812106687651020584; 364 binterpt[1][1]=0.559074397713145440020984353917; 365 binterpt[2][1]=0.16666666666666666666666666666667; 366 367 ierr = TSRosWRegister(TSROSWSANDU3,3,3,&A[0][0],&Gamma[0][0],b,b2,2,&binterpt[0][0]);CHKERRQ(ierr); 368 } 369 { 370 const PetscReal s3 = PetscSqrtReal(3.),g = (3.0+s3)/6.0; 371 const PetscReal 372 A[3][3] = {{0,0,0}, 373 {1,0,0}, 374 {0.25,0.25,0}}, 375 Gamma[3][3] = {{0,0,0}, 376 {(-3.0-s3)/6.0,g,0}, 377 {(-3.0-s3)/24.0,(-3.0-s3)/8.0,g}}, 378 b[3] = {1./6.,1./6.,2./3.}, 379 b2[3] = {1./4.,1./4.,1./2.}; 380 381 PetscReal binterpt[3][2]; 382 binterpt[0][0]=0.089316397477040902157517886164709; 383 binterpt[1][0]=-0.91068360252295909784248211383529; 384 binterpt[2][0]=1.8213672050459181956849642276706; 385 binterpt[0][1]=0.077350269189625764509148780501957; 386 binterpt[1][1]=1.077350269189625764509148780502; 387 binterpt[2][1]=-1.1547005383792515290182975610039; 388 ierr = TSRosWRegister(TSROSWASSP3P3S1C,3,3,&A[0][0],&Gamma[0][0],b,b2,2,&binterpt[0][0]);CHKERRQ(ierr); 389 } 390 391 { 392 const PetscReal 393 A[4][4] = {{0,0,0,0}, 394 {1./2.,0,0,0}, 395 {1./2.,1./2.,0,0}, 396 {1./6.,1./6.,1./6.,0}}, 397 Gamma[4][4] = {{1./2.,0,0,0}, 398 {0.0,1./4.,0,0}, 399 {-2.,-2./3.,2./3.,0}, 400 {1./2.,5./36.,-2./9,0}}, 401 b[4] = {1./6.,1./6.,1./6.,1./2.}, 402 b2[4] = {1./8.,3./4.,1./8.,0}; 403 PetscReal binterpt[4][3]; 404 binterpt[0][0]=6.25; 405 binterpt[1][0]=-30.25; 406 binterpt[2][0]=1.75; 407 binterpt[3][0]=23.25; 408 binterpt[0][1]=-9.75; 409 binterpt[1][1]=58.75; 410 binterpt[2][1]=-3.25; 411 binterpt[3][1]=-45.75; 412 binterpt[0][2]=3.6666666666666666666666666666667; 413 binterpt[1][2]=-28.333333333333333333333333333333; 414 binterpt[2][2]=1.6666666666666666666666666666667; 415 binterpt[3][2]=23.; 416 ierr = TSRosWRegister(TSROSWLASSP3P4S2C,3,4,&A[0][0],&Gamma[0][0],b,b2,3,&binterpt[0][0]);CHKERRQ(ierr); 417 } 418 419 { 420 const PetscReal 421 A[4][4] = {{0,0,0,0}, 422 {1./2.,0,0,0}, 423 {1./2.,1./2.,0,0}, 424 {1./6.,1./6.,1./6.,0}}, 425 Gamma[4][4] = {{1./2.,0,0,0}, 426 {0.0,3./4.,0,0}, 427 {-2./3.,-23./9.,2./9.,0}, 428 {1./18.,65./108.,-2./27,0}}, 429 b[4] = {1./6.,1./6.,1./6.,1./2.}, 430 b2[4] = {3./16.,10./16.,3./16.,0}; 431 432 PetscReal binterpt[4][3]; 433 binterpt[0][0]=1.6911764705882352941176470588235; 434 binterpt[1][0]=3.6813725490196078431372549019608; 435 binterpt[2][0]=0.23039215686274509803921568627451; 436 binterpt[3][0]=-4.6029411764705882352941176470588; 437 binterpt[0][1]=-0.95588235294117647058823529411765; 438 binterpt[1][1]=-6.2401960784313725490196078431373; 439 binterpt[2][1]=-0.31862745098039215686274509803922; 440 binterpt[3][1]=7.5147058823529411764705882352941; 441 binterpt[0][2]=-0.56862745098039215686274509803922; 442 binterpt[1][2]=2.7254901960784313725490196078431; 443 binterpt[2][2]=0.25490196078431372549019607843137; 444 binterpt[3][2]=-2.4117647058823529411764705882353; 445 ierr = TSRosWRegister(TSROSWLLSSP3P4S2C,3,4,&A[0][0],&Gamma[0][0],b,b2,3,&binterpt[0][0]);CHKERRQ(ierr); 446 } 447 448 { 449 PetscReal A[4][4],Gamma[4][4],b[4],b2[4]; 450 PetscReal binterpt[4][3]; 451 452 Gamma[0][0]=0.4358665215084589994160194475295062513822671686978816; 453 Gamma[0][1]=0; Gamma[0][2]=0; Gamma[0][3]=0; 454 Gamma[1][0]=-1.997527830934941248426324674704153457289527280554476; 455 Gamma[1][1]=0.4358665215084589994160194475295062513822671686978816; 456 Gamma[1][2]=0; Gamma[1][3]=0; 457 Gamma[2][0]=-1.007948511795029620852002345345404191008352770119903; 458 Gamma[2][1]=-0.004648958462629345562774289390054679806993396798458131; 459 Gamma[2][2]=0.4358665215084589994160194475295062513822671686978816; 460 Gamma[2][3]=0; 461 Gamma[3][0]=-0.6685429734233467180451604600279552604364311322650783; 462 Gamma[3][1]=0.6056625986449338476089525334450053439525178740492984; 463 Gamma[3][2]=-0.9717899277217721234705114616271378792182450260943198; 464 Gamma[3][3]=0; 465 466 A[0][0]=0; A[0][1]=0; A[0][2]=0; A[0][3]=0; 467 A[1][0]=0.8717330430169179988320388950590125027645343373957631; 468 A[1][1]=0; A[1][2]=0; A[1][3]=0; 469 A[2][0]=0.5275890119763004115618079766722914408876108660811028; 470 A[2][1]=0.07241098802369958843819203208518599088698057726988732; 471 A[2][2]=0; A[2][3]=0; 472 A[3][0]=0.3990960076760701320627260685975778145384666450351314; 473 A[3][1]=-0.4375576546135194437228463747348862825846903771419953; 474 A[3][2]=1.038461646937449311660120300601880176655352737312713; 475 A[3][3]=0; 476 477 b[0]=0.1876410243467238251612921333138006734899663569186926; 478 b[1]=-0.5952974735769549480478230473706443582188442040780541; 479 b[2]=0.9717899277217721234705114616271378792182450260943198; 480 b[3]=0.4358665215084589994160194475295062513822671686978816; 481 482 b2[0]=0.2147402862233891404862383521089097657790734483804460; 483 b2[1]=-0.4851622638849390928209050538171743017757490232519684; 484 b2[2]=0.8687250025203875511662123688667549217531982787600080; 485 b2[3]=0.4016969751411624011684543450940068201770721128357014; 486 487 binterpt[0][0]=2.2565812720167954547104627844105; 488 binterpt[1][0]=1.349166413351089573796243820819; 489 binterpt[2][0]=-2.4695174540533503758652847586647; 490 binterpt[3][0]=-0.13623023131453465264142184656474; 491 binterpt[0][1]=-3.0826699111559187902922463354557; 492 binterpt[1][1]=-2.4689115685996042534544925650515; 493 binterpt[2][1]=5.7428279814696677152129332773553; 494 binterpt[3][1]=-0.19124650171414467146619437684812; 495 binterpt[0][2]=1.0137296634858471607430756831148; 496 binterpt[1][2]=0.52444768167155973161042570784064; 497 binterpt[2][2]=-2.3015205996945452158771370439586; 498 binterpt[3][2]=0.76334325453713832352363565300308; 499 500 ierr = TSRosWRegister(TSROSWARK3,3,4,&A[0][0],&Gamma[0][0],b,b2,3,&binterpt[0][0]);CHKERRQ(ierr); 501 } 502 503 PetscFunctionReturn(0); 504 } 505 506 #undef __FUNCT__ 507 #define __FUNCT__ "TSRosWRegisterDestroy" 508 /*@C 509 TSRosWRegisterDestroy - Frees the list of schemes that were registered by TSRosWRegister(). 510 511 Not Collective 512 513 Level: advanced 514 515 .keywords: TSRosW, register, destroy 516 .seealso: TSRosWRegister(), TSRosWRegisterAll(), TSRosWRegisterDynamic() 517 @*/ 518 PetscErrorCode TSRosWRegisterDestroy(void) 519 { 520 PetscErrorCode ierr; 521 RosWTableauLink link; 522 523 PetscFunctionBegin; 524 while ((link = RosWTableauList)) { 525 RosWTableau t = &link->tab; 526 RosWTableauList = link->next; 527 ierr = PetscFree5(t->A,t->Gamma,t->b,t->ASum,t->GammaSum);CHKERRQ(ierr); 528 ierr = PetscFree5(t->At,t->bt,t->GammaInv,t->GammaZeroDiag,t->GammaExplicitCorr);CHKERRQ(ierr); 529 ierr = PetscFree2(t->bembed,t->bembedt);CHKERRQ(ierr); 530 ierr = PetscFree(t->binterpt);CHKERRQ(ierr); 531 ierr = PetscFree(t->name);CHKERRQ(ierr); 532 ierr = PetscFree(link);CHKERRQ(ierr); 533 } 534 TSRosWRegisterAllCalled = PETSC_FALSE; 535 PetscFunctionReturn(0); 536 } 537 538 #undef __FUNCT__ 539 #define __FUNCT__ "TSRosWInitializePackage" 540 /*@C 541 TSRosWInitializePackage - This function initializes everything in the TSRosW package. It is called 542 from PetscDLLibraryRegister() when using dynamic libraries, and on the first call to TSCreate_RosW() 543 when using static libraries. 544 545 Input Parameter: 546 path - The dynamic library path, or PETSC_NULL 547 548 Level: developer 549 550 .keywords: TS, TSRosW, initialize, package 551 .seealso: PetscInitialize() 552 @*/ 553 PetscErrorCode TSRosWInitializePackage(const char path[]) 554 { 555 PetscErrorCode ierr; 556 557 PetscFunctionBegin; 558 if (TSRosWPackageInitialized) PetscFunctionReturn(0); 559 TSRosWPackageInitialized = PETSC_TRUE; 560 ierr = TSRosWRegisterAll();CHKERRQ(ierr); 561 ierr = PetscRegisterFinalize(TSRosWFinalizePackage);CHKERRQ(ierr); 562 PetscFunctionReturn(0); 563 } 564 565 #undef __FUNCT__ 566 #define __FUNCT__ "TSRosWFinalizePackage" 567 /*@C 568 TSRosWFinalizePackage - This function destroys everything in the TSRosW package. It is 569 called from PetscFinalize(). 570 571 Level: developer 572 573 .keywords: Petsc, destroy, package 574 .seealso: PetscFinalize() 575 @*/ 576 PetscErrorCode TSRosWFinalizePackage(void) 577 { 578 PetscErrorCode ierr; 579 580 PetscFunctionBegin; 581 TSRosWPackageInitialized = PETSC_FALSE; 582 ierr = TSRosWRegisterDestroy();CHKERRQ(ierr); 583 PetscFunctionReturn(0); 584 } 585 586 #undef __FUNCT__ 587 #define __FUNCT__ "TSRosWRegister" 588 /*@C 589 TSRosWRegister - register a Rosenbrock W scheme by providing the entries in the Butcher tableau and optionally embedded approximations and interpolation 590 591 Not Collective, but the same schemes should be registered on all processes on which they will be used 592 593 Input Parameters: 594 + name - identifier for method 595 . order - approximation order of method 596 . s - number of stages, this is the dimension of the matrices below 597 . A - Table of propagated stage coefficients (dimension s*s, row-major), strictly lower triangular 598 . Gamma - Table of coefficients in implicit stage equations (dimension s*s, row-major), lower triangular with nonzero diagonal 599 . b - Step completion table (dimension s) 600 - bembed - Step completion table for a scheme of order one less (dimension s, PETSC_NULL if no embedded scheme is available) 601 . pinterp - Order of the interpolation scheme, equal to the number of columns of binterpt 602 . binterpt - Coefficients of the interpolation formula (dimension s*pinterp) 603 604 Notes: 605 Several Rosenbrock W methods are provided, this function is only needed to create new methods. 606 607 Level: advanced 608 609 .keywords: TS, register 610 611 .seealso: TSRosW 612 @*/ 613 PetscErrorCode TSRosWRegister(const TSRosWType name,PetscInt order,PetscInt s, 614 const PetscReal A[],const PetscReal Gamma[],const PetscReal b[],const PetscReal bembed[], 615 PetscInt pinterp,const PetscReal binterpt[]) 616 { 617 PetscErrorCode ierr; 618 RosWTableauLink link; 619 RosWTableau t; 620 PetscInt i,j,k; 621 PetscScalar *GammaInv; 622 623 PetscFunctionBegin; 624 PetscValidCharPointer(name,1); 625 PetscValidPointer(A,4); 626 PetscValidPointer(Gamma,5); 627 PetscValidPointer(b,6); 628 if (bembed) PetscValidPointer(bembed,7); 629 630 ierr = PetscMalloc(sizeof(*link),&link);CHKERRQ(ierr); 631 ierr = PetscMemzero(link,sizeof(*link));CHKERRQ(ierr); 632 t = &link->tab; 633 ierr = PetscStrallocpy(name,&t->name);CHKERRQ(ierr); 634 t->order = order; 635 t->s = s; 636 ierr = PetscMalloc5(s*s,PetscReal,&t->A,s*s,PetscReal,&t->Gamma,s,PetscReal,&t->b,s,PetscReal,&t->ASum,s,PetscReal,&t->GammaSum);CHKERRQ(ierr); 637 ierr = PetscMalloc5(s*s,PetscReal,&t->At,s,PetscReal,&t->bt,s*s,PetscReal,&t->GammaInv,s,PetscBool,&t->GammaZeroDiag,s*s,PetscReal,&t->GammaExplicitCorr);CHKERRQ(ierr); 638 ierr = PetscMemcpy(t->A,A,s*s*sizeof(A[0]));CHKERRQ(ierr); 639 ierr = PetscMemcpy(t->Gamma,Gamma,s*s*sizeof(Gamma[0]));CHKERRQ(ierr); 640 ierr = PetscMemcpy(t->GammaExplicitCorr,Gamma,s*s*sizeof(Gamma[0]));CHKERRQ(ierr); 641 ierr = PetscMemcpy(t->b,b,s*sizeof(b[0]));CHKERRQ(ierr); 642 if (bembed) { 643 ierr = PetscMalloc2(s,PetscReal,&t->bembed,s,PetscReal,&t->bembedt);CHKERRQ(ierr); 644 ierr = PetscMemcpy(t->bembed,bembed,s*sizeof(bembed[0]));CHKERRQ(ierr); 645 } 646 for (i=0; i<s; i++) { 647 t->ASum[i] = 0; 648 t->GammaSum[i] = 0; 649 for (j=0; j<s; j++) { 650 t->ASum[i] += A[i*s+j]; 651 t->GammaSum[i] += Gamma[i*s+j]; 652 } 653 } 654 ierr = PetscMalloc(s*s*sizeof(PetscScalar),&GammaInv);CHKERRQ(ierr); /* Need to use Scalar for inverse, then convert back to Real */ 655 for (i=0; i<s*s; i++) GammaInv[i] = Gamma[i]; 656 for (i=0; i<s; i++) { 657 if (Gamma[i*s+i] == 0.0) { 658 GammaInv[i*s+i] = 1.0; 659 t->GammaZeroDiag[i] = PETSC_TRUE; 660 } else { 661 t->GammaZeroDiag[i] = PETSC_FALSE; 662 } 663 } 664 665 switch (s) { 666 case 1: GammaInv[0] = 1./GammaInv[0]; break; 667 case 2: ierr = PetscKernel_A_gets_inverse_A_2(GammaInv,0);CHKERRQ(ierr); break; 668 case 3: ierr = PetscKernel_A_gets_inverse_A_3(GammaInv,0);CHKERRQ(ierr); break; 669 case 4: ierr = PetscKernel_A_gets_inverse_A_4(GammaInv,0);CHKERRQ(ierr); break; 670 case 5: { 671 PetscInt ipvt5[5]; 672 MatScalar work5[5*5]; 673 ierr = PetscKernel_A_gets_inverse_A_5(GammaInv,ipvt5,work5,0);CHKERRQ(ierr); break; 674 } 675 case 6: ierr = PetscKernel_A_gets_inverse_A_6(GammaInv,0);CHKERRQ(ierr); break; 676 case 7: ierr = PetscKernel_A_gets_inverse_A_7(GammaInv,0);CHKERRQ(ierr); break; 677 default: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not implemented for %D stages",s); 678 } 679 for (i=0; i<s*s; i++) t->GammaInv[i] = PetscRealPart(GammaInv[i]); 680 ierr = PetscFree(GammaInv);CHKERRQ(ierr); 681 682 for (i=0; i<s; i++) { 683 for (k=0; k<i+1; k++) { 684 t->GammaExplicitCorr[i*s+k]=(t->GammaExplicitCorr[i*s+k])*(t->GammaInv[k*s+k]); 685 for (j=k+1; j<i+1; j++) { 686 t->GammaExplicitCorr[i*s+k]+=(t->GammaExplicitCorr[i*s+j])*(t->GammaInv[j*s+k]); 687 } 688 } 689 } 690 691 for (i=0; i<s; i++) { 692 for (j=0; j<s; j++) { 693 t->At[i*s+j] = 0; 694 for (k=0; k<s; k++) { 695 t->At[i*s+j] += t->A[i*s+k] * t->GammaInv[k*s+j]; 696 } 697 } 698 t->bt[i] = 0; 699 for (j=0; j<s; j++) { 700 t->bt[i] += t->b[j] * t->GammaInv[j*s+i]; 701 } 702 if (bembed) { 703 t->bembedt[i] = 0; 704 for (j=0; j<s; j++) { 705 t->bembedt[i] += t->bembed[j] * t->GammaInv[j*s+i]; 706 } 707 } 708 } 709 t->ccfl = 1.0; /* Fix this */ 710 711 t->pinterp = pinterp; 712 ierr = PetscMalloc(s*pinterp*sizeof(binterpt[0]),&t->binterpt);CHKERRQ(ierr); 713 ierr = PetscMemcpy(t->binterpt,binterpt,s*pinterp*sizeof(binterpt[0]));CHKERRQ(ierr); 714 link->next = RosWTableauList; 715 RosWTableauList = link; 716 PetscFunctionReturn(0); 717 } 718 719 #undef __FUNCT__ 720 #define __FUNCT__ "TSEvaluateStep_RosW" 721 /* 722 The step completion formula is 723 724 x1 = x0 + b^T Y 725 726 where Y is the multi-vector of stages corrections. This function can be called before or after ts->vec_sol has been 727 updated. Suppose we have a completion formula b and an embedded formula be of different order. We can write 728 729 x1e = x0 + be^T Y 730 = x1 - b^T Y + be^T Y 731 = x1 + (be - b)^T Y 732 733 so we can evaluate the method of different order even after the step has been optimistically completed. 734 */ 735 static PetscErrorCode TSEvaluateStep_RosW(TS ts,PetscInt order,Vec X,PetscBool *done) 736 { 737 TS_RosW *ros = (TS_RosW*)ts->data; 738 RosWTableau tab = ros->tableau; 739 PetscScalar *w = ros->work; 740 PetscInt i; 741 PetscErrorCode ierr; 742 743 PetscFunctionBegin; 744 if (order == tab->order) { 745 if (ros->status == TS_STEP_INCOMPLETE) { /* Use standard completion formula */ 746 ierr = VecCopy(ts->vec_sol,X);CHKERRQ(ierr); 747 for (i=0; i<tab->s; i++) w[i] = tab->bt[i]; 748 ierr = VecMAXPY(X,tab->s,w,ros->Y);CHKERRQ(ierr); 749 } else {ierr = VecCopy(ts->vec_sol,X);CHKERRQ(ierr);} 750 if (done) *done = PETSC_TRUE; 751 PetscFunctionReturn(0); 752 } else if (order == tab->order-1) { 753 if (!tab->bembedt) goto unavailable; 754 if (ros->status == TS_STEP_INCOMPLETE) { /* Use embedded completion formula */ 755 ierr = VecCopy(ts->vec_sol,X);CHKERRQ(ierr); 756 for (i=0; i<tab->s; i++) w[i] = tab->bembedt[i]; 757 ierr = VecMAXPY(X,tab->s,w,ros->Y);CHKERRQ(ierr); 758 } else { /* Use rollback-and-recomplete formula (bembedt - bt) */ 759 for (i=0; i<tab->s; i++) w[i] = tab->bembedt[i] - tab->bt[i]; 760 ierr = VecCopy(ts->vec_sol,X);CHKERRQ(ierr); 761 ierr = VecMAXPY(X,tab->s,w,ros->Y);CHKERRQ(ierr); 762 } 763 if (done) *done = PETSC_TRUE; 764 PetscFunctionReturn(0); 765 } 766 unavailable: 767 if (done) *done = PETSC_FALSE; 768 else SETERRQ3(((PetscObject)ts)->comm,PETSC_ERR_SUP,"Rosenbrock-W '%s' of order %D cannot evaluate step at order %D",tab->name,tab->order,order); 769 PetscFunctionReturn(0); 770 } 771 772 #undef __FUNCT__ 773 #define __FUNCT__ "TSStep_RosW" 774 static PetscErrorCode TSStep_RosW(TS ts) 775 { 776 TS_RosW *ros = (TS_RosW*)ts->data; 777 RosWTableau tab = ros->tableau; 778 const PetscInt s = tab->s; 779 const PetscReal *At = tab->At,*Gamma = tab->Gamma,*ASum = tab->ASum,*GammaInv = tab->GammaInv; 780 const PetscReal *GammaExplicitCorr = tab->GammaExplicitCorr; 781 const PetscBool *GammaZeroDiag = tab->GammaZeroDiag; 782 PetscScalar *w = ros->work; 783 Vec *Y = ros->Y,Ydot = ros->Ydot,Zdot = ros->Zdot,Zstage = ros->Zstage; 784 SNES snes; 785 TSAdapt adapt; 786 PetscInt i,j,its,lits,reject,next_scheme; 787 PetscReal next_time_step; 788 PetscBool accept; 789 PetscErrorCode ierr; 790 MatStructure str; 791 792 PetscFunctionBegin; 793 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 794 next_time_step = ts->time_step; 795 accept = PETSC_TRUE; 796 ros->status = TS_STEP_INCOMPLETE; 797 798 for (reject=0; reject<ts->max_reject && !ts->reason; reject++,ts->reject++) { 799 const PetscReal h = ts->time_step; 800 ierr = TSPreStep(ts);CHKERRQ(ierr); 801 ierr = VecCopy(ts->vec_sol,ros->VecSolPrev);CHKERRQ(ierr);/*move this at the end*/ 802 for (i=0; i<s; i++) { 803 ros->stage_time = ts->ptime + h*ASum[i]; 804 ierr = TSPreStage(ts,ros->stage_time);CHKERRQ(ierr); 805 if (GammaZeroDiag[i]) { 806 ros->stage_explicit = PETSC_TRUE; 807 ros->shift = 1./h; 808 } else { 809 ros->stage_explicit = PETSC_FALSE; 810 ros->shift = 1./(h*Gamma[i*s+i]); 811 } 812 813 ierr = VecCopy(ts->vec_sol,Zstage);CHKERRQ(ierr); 814 for (j=0; j<i; j++) w[j] = At[i*s+j]; 815 ierr = VecMAXPY(Zstage,i,w,Y);CHKERRQ(ierr); 816 817 for (j=0; j<i; j++) w[j] = 1./h * GammaInv[i*s+j]; 818 ierr = VecZeroEntries(Zdot);CHKERRQ(ierr); 819 ierr = VecMAXPY(Zdot,i,w,Y);CHKERRQ(ierr); 820 821 /* Initial guess taken from last stage */ 822 ierr = VecZeroEntries(Y[i]);CHKERRQ(ierr); 823 824 if (!ros->stage_explicit) { 825 if (!ros->recompute_jacobian && !i) { 826 ierr = SNESSetLagJacobian(snes,-2);CHKERRQ(ierr); /* Recompute the Jacobian on this solve, but not again */ 827 } 828 ierr = SNESSolve(snes,PETSC_NULL,Y[i]);CHKERRQ(ierr); 829 ierr = SNESGetIterationNumber(snes,&its);CHKERRQ(ierr); 830 ierr = SNESGetLinearSolveIterations(snes,&lits);CHKERRQ(ierr); 831 ts->snes_its += its; ts->ksp_its += lits; 832 ierr = TSGetAdapt(ts,&adapt);CHKERRQ(ierr); 833 ierr = TSAdaptCheckStage(adapt,ts,&accept);CHKERRQ(ierr); 834 if (!accept) goto reject_step; 835 } else { 836 Mat J,Jp; 837 ierr = VecZeroEntries(Ydot);CHKERRQ(ierr); /* Evaluate Y[i]=G(t,Ydot=0,Zstage) */ 838 ierr = TSComputeIFunction(ts,ros->stage_time,Zstage,Ydot,Y[i],PETSC_FALSE);CHKERRQ(ierr); 839 ierr = VecScale(Y[i],-1.0); 840 ierr = VecAXPY(Y[i],-1.0,Zdot);CHKERRQ(ierr); /*Y[i]=F(Zstage)-Zdot[=GammaInv*Y]*/ 841 842 ierr = VecZeroEntries(Zstage);CHKERRQ(ierr); /* Zstage = GammaExplicitCorr[i,j] * Y[j] */ 843 for (j=0; j<i; j++) w[j] = GammaExplicitCorr[i*s+j]; 844 ierr = VecMAXPY(Zstage,i,w,Y);CHKERRQ(ierr); 845 /*Y[i] += Y[i] + Jac*Zstage[=Jac*GammaExplicitCorr[i,j] * Y[j]] */ 846 str = SAME_NONZERO_PATTERN; 847 ierr = TSGetIJacobian(ts,&J,&Jp,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr); 848 ierr = TSComputeIJacobian(ts,ros->stage_time,ts->vec_sol,Ydot,0,&J,&Jp,&str,PETSC_FALSE);CHKERRQ(ierr); 849 ierr = MatMult(J,Zstage,Zdot); 850 851 ierr = VecAXPY(Y[i],-1.0,Zdot);CHKERRQ(ierr); 852 ierr = VecScale(Y[i],h); 853 ts->ksp_its += 1; 854 } 855 } 856 ierr = TSEvaluateStep(ts,tab->order,ts->vec_sol,PETSC_NULL);CHKERRQ(ierr); 857 ros->status = TS_STEP_PENDING; 858 859 /* Register only the current method as a candidate because we're not supporting multiple candidates yet. */ 860 ierr = TSGetAdapt(ts,&adapt);CHKERRQ(ierr); 861 ierr = TSAdaptCandidatesClear(adapt);CHKERRQ(ierr); 862 ierr = TSAdaptCandidateAdd(adapt,tab->name,tab->order,1,tab->ccfl,1.*tab->s,PETSC_TRUE);CHKERRQ(ierr); 863 ierr = TSAdaptChoose(adapt,ts,ts->time_step,&next_scheme,&next_time_step,&accept);CHKERRQ(ierr); 864 if (accept) { 865 /* ignore next_scheme for now */ 866 ts->ptime += ts->time_step; 867 ts->time_step = next_time_step; 868 ts->steps++; 869 ros->status = TS_STEP_COMPLETE; 870 break; 871 } else { /* Roll back the current step */ 872 for (i=0; i<s; i++) w[i] = -tab->bt[i]; 873 ierr = VecMAXPY(ts->vec_sol,s,w,Y);CHKERRQ(ierr); 874 ts->time_step = next_time_step; 875 ros->status = TS_STEP_INCOMPLETE; 876 } 877 reject_step: continue; 878 } 879 if (ros->status != TS_STEP_COMPLETE && !ts->reason) ts->reason = TS_DIVERGED_STEP_REJECTED; 880 PetscFunctionReturn(0); 881 } 882 883 #undef __FUNCT__ 884 #define __FUNCT__ "TSInterpolate_RosW" 885 static PetscErrorCode TSInterpolate_RosW(TS ts,PetscReal itime,Vec X) 886 { 887 TS_RosW *ros = (TS_RosW*)ts->data; 888 PetscInt s = ros->tableau->s,pinterp = ros->tableau->pinterp,i,j; 889 PetscReal h; 890 PetscReal tt,t; 891 PetscScalar *bt; 892 const PetscReal *Bt = ros->tableau->binterpt; 893 PetscErrorCode ierr; 894 const PetscReal *GammaInv = ros->tableau->GammaInv; 895 PetscScalar *w = ros->work; 896 Vec *Y = ros->Y; 897 898 PetscFunctionBegin; 899 if (!Bt) SETERRQ1(((PetscObject)ts)->comm,PETSC_ERR_SUP,"TSRosW %s does not have an interpolation formula",ros->tableau->name); 900 901 switch (ros->status) { 902 case TS_STEP_INCOMPLETE: 903 case TS_STEP_PENDING: 904 h = ts->time_step; 905 t = (itime - ts->ptime)/h; 906 break; 907 case TS_STEP_COMPLETE: 908 h = ts->time_step_prev; 909 t = (itime - ts->ptime)/h + 1; /* In the interval [0,1] */ 910 break; 911 default: SETERRQ(((PetscObject)ts)->comm,PETSC_ERR_PLIB,"Invalid TSStepStatus"); 912 } 913 ierr = PetscMalloc(s*sizeof(bt[0]),&bt);CHKERRQ(ierr); 914 for (i=0; i<s; i++) bt[i] = 0; 915 for (j=0,tt=t; j<pinterp; j++,tt*=t) { 916 for (i=0; i<s; i++) { 917 bt[i] += Bt[i*pinterp+j] * tt; 918 } 919 } 920 921 /* y(t+tt*h) = y(t) + Sum bt(tt) * GammaInv * Ydot */ 922 /*X<-0*/ 923 ierr = VecZeroEntries(X);CHKERRQ(ierr); 924 925 /*X<- Sum bt_i * GammaInv(i,1:i) * Y(1:i) */ 926 for (j=0; j<s; j++) w[j]=0; 927 for (j=0; j<s; j++) { 928 for (i=j; i<s; i++) { 929 w[j] += bt[i]*GammaInv[i*s+j]; 930 } 931 } 932 ierr = VecMAXPY(X,i,w,Y);CHKERRQ(ierr); 933 934 /*X<-y(t) + X*/ 935 ierr = VecAXPY(X,1.0,ros->VecSolPrev);CHKERRQ(ierr); 936 937 ierr = PetscFree(bt);CHKERRQ(ierr); 938 939 PetscFunctionReturn(0); 940 } 941 942 /*------------------------------------------------------------*/ 943 #undef __FUNCT__ 944 #define __FUNCT__ "TSReset_RosW" 945 static PetscErrorCode TSReset_RosW(TS ts) 946 { 947 TS_RosW *ros = (TS_RosW*)ts->data; 948 PetscInt s; 949 PetscErrorCode ierr; 950 951 PetscFunctionBegin; 952 if (!ros->tableau) PetscFunctionReturn(0); 953 s = ros->tableau->s; 954 ierr = VecDestroyVecs(s,&ros->Y);CHKERRQ(ierr); 955 ierr = VecDestroy(&ros->Ydot);CHKERRQ(ierr); 956 ierr = VecDestroy(&ros->Ystage);CHKERRQ(ierr); 957 ierr = VecDestroy(&ros->Zdot);CHKERRQ(ierr); 958 ierr = VecDestroy(&ros->Zstage);CHKERRQ(ierr); 959 ierr = VecDestroy(&ros->VecSolPrev);CHKERRQ(ierr); 960 ierr = PetscFree(ros->work);CHKERRQ(ierr); 961 PetscFunctionReturn(0); 962 } 963 964 #undef __FUNCT__ 965 #define __FUNCT__ "TSDestroy_RosW" 966 static PetscErrorCode TSDestroy_RosW(TS ts) 967 { 968 PetscErrorCode ierr; 969 970 PetscFunctionBegin; 971 ierr = TSReset_RosW(ts);CHKERRQ(ierr); 972 ierr = PetscFree(ts->data);CHKERRQ(ierr); 973 ierr = PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSRosWGetType_C","",PETSC_NULL);CHKERRQ(ierr); 974 ierr = PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSRosWSetType_C","",PETSC_NULL);CHKERRQ(ierr); 975 ierr = PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSRosWSetRecomputeJacobian_C","",PETSC_NULL);CHKERRQ(ierr); 976 PetscFunctionReturn(0); 977 } 978 979 /* 980 This defines the nonlinear equation that is to be solved with SNES 981 G(U) = F[t0+Theta*dt, U, (U-U0)*shift] = 0 982 */ 983 #undef __FUNCT__ 984 #define __FUNCT__ "SNESTSFormFunction_RosW" 985 static PetscErrorCode SNESTSFormFunction_RosW(SNES snes,Vec X,Vec F,TS ts) 986 { 987 TS_RosW *ros = (TS_RosW*)ts->data; 988 PetscErrorCode ierr; 989 990 PetscFunctionBegin; 991 ierr = VecWAXPY(ros->Ydot,ros->shift,X,ros->Zdot);CHKERRQ(ierr); /* Ydot = shift*X + Zdot */ 992 ierr = VecWAXPY(ros->Ystage,1.0,X,ros->Zstage);CHKERRQ(ierr); /* Ystage = X + Zstage */ 993 ierr = TSComputeIFunction(ts,ros->stage_time,ros->Ystage,ros->Ydot,F,PETSC_FALSE);CHKERRQ(ierr); 994 PetscFunctionReturn(0); 995 } 996 997 #undef __FUNCT__ 998 #define __FUNCT__ "SNESTSFormJacobian_RosW" 999 static PetscErrorCode SNESTSFormJacobian_RosW(SNES snes,Vec X,Mat *A,Mat *B,MatStructure *str,TS ts) 1000 { 1001 TS_RosW *ros = (TS_RosW*)ts->data; 1002 PetscErrorCode ierr; 1003 1004 PetscFunctionBegin; 1005 /* ros->Ydot and ros->Ystage have already been computed in SNESTSFormFunction_RosW (SNES guarantees this) */ 1006 ierr = TSComputeIJacobian(ts,ros->stage_time,ros->Ystage,ros->Ydot,ros->shift,A,B,str,PETSC_TRUE);CHKERRQ(ierr); 1007 PetscFunctionReturn(0); 1008 } 1009 1010 #undef __FUNCT__ 1011 #define __FUNCT__ "TSSetUp_RosW" 1012 static PetscErrorCode TSSetUp_RosW(TS ts) 1013 { 1014 TS_RosW *ros = (TS_RosW*)ts->data; 1015 RosWTableau tab = ros->tableau; 1016 PetscInt s = tab->s; 1017 PetscErrorCode ierr; 1018 1019 PetscFunctionBegin; 1020 if (!ros->tableau) { 1021 ierr = TSRosWSetType(ts,TSRosWDefault);CHKERRQ(ierr); 1022 } 1023 ierr = VecDuplicateVecs(ts->vec_sol,s,&ros->Y);CHKERRQ(ierr); 1024 ierr = VecDuplicate(ts->vec_sol,&ros->Ydot);CHKERRQ(ierr); 1025 ierr = VecDuplicate(ts->vec_sol,&ros->Ystage);CHKERRQ(ierr); 1026 ierr = VecDuplicate(ts->vec_sol,&ros->Zdot);CHKERRQ(ierr); 1027 ierr = VecDuplicate(ts->vec_sol,&ros->Zstage);CHKERRQ(ierr); 1028 ierr = VecDuplicate(ts->vec_sol,&ros->VecSolPrev);CHKERRQ(ierr); 1029 ierr = PetscMalloc(s*sizeof(ros->work[0]),&ros->work);CHKERRQ(ierr); 1030 PetscFunctionReturn(0); 1031 } 1032 /*------------------------------------------------------------*/ 1033 1034 #undef __FUNCT__ 1035 #define __FUNCT__ "TSSetFromOptions_RosW" 1036 static PetscErrorCode TSSetFromOptions_RosW(TS ts) 1037 { 1038 TS_RosW *ros = (TS_RosW*)ts->data; 1039 PetscErrorCode ierr; 1040 char rostype[256]; 1041 1042 PetscFunctionBegin; 1043 ierr = PetscOptionsHead("RosW ODE solver options");CHKERRQ(ierr); 1044 { 1045 RosWTableauLink link; 1046 PetscInt count,choice; 1047 PetscBool flg; 1048 const char **namelist; 1049 SNES snes; 1050 1051 ierr = PetscStrncpy(rostype,TSRosWDefault,sizeof rostype);CHKERRQ(ierr); 1052 for (link=RosWTableauList,count=0; link; link=link->next,count++) ; 1053 ierr = PetscMalloc(count*sizeof(char*),&namelist);CHKERRQ(ierr); 1054 for (link=RosWTableauList,count=0; link; link=link->next,count++) namelist[count] = link->tab.name; 1055 ierr = PetscOptionsEList("-ts_rosw_type","Family of Rosenbrock-W method","TSRosWSetType",(const char*const*)namelist,count,rostype,&choice,&flg);CHKERRQ(ierr); 1056 ierr = TSRosWSetType(ts,flg ? namelist[choice] : rostype);CHKERRQ(ierr); 1057 ierr = PetscFree(namelist);CHKERRQ(ierr); 1058 1059 ierr = PetscOptionsBool("-ts_rosw_recompute_jacobian","Recompute the Jacobian at each stage","TSRosWSetRecomputeJacobian",ros->recompute_jacobian,&ros->recompute_jacobian,PETSC_NULL);CHKERRQ(ierr); 1060 1061 /* Rosenbrock methods are linearly implicit, so set that unless the user has specifically asked for something else */ 1062 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 1063 if (!((PetscObject)snes)->type_name) { 1064 ierr = SNESSetType(snes,SNESKSPONLY);CHKERRQ(ierr); 1065 } 1066 ierr = SNESSetFromOptions(snes);CHKERRQ(ierr); 1067 } 1068 ierr = PetscOptionsTail();CHKERRQ(ierr); 1069 PetscFunctionReturn(0); 1070 } 1071 1072 #undef __FUNCT__ 1073 #define __FUNCT__ "PetscFormatRealArray" 1074 static PetscErrorCode PetscFormatRealArray(char buf[],size_t len,const char *fmt,PetscInt n,const PetscReal x[]) 1075 { 1076 PetscErrorCode ierr; 1077 PetscInt i; 1078 size_t left,count; 1079 char *p; 1080 1081 PetscFunctionBegin; 1082 for (i=0,p=buf,left=len; i<n; i++) { 1083 ierr = PetscSNPrintfCount(p,left,fmt,&count,x[i]);CHKERRQ(ierr); 1084 if (count >= left) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Insufficient space in buffer"); 1085 left -= count; 1086 p += count; 1087 *p++ = ' '; 1088 } 1089 p[i ? 0 : -1] = 0; 1090 PetscFunctionReturn(0); 1091 } 1092 1093 #undef __FUNCT__ 1094 #define __FUNCT__ "TSView_RosW" 1095 static PetscErrorCode TSView_RosW(TS ts,PetscViewer viewer) 1096 { 1097 TS_RosW *ros = (TS_RosW*)ts->data; 1098 RosWTableau tab = ros->tableau; 1099 PetscBool iascii; 1100 PetscErrorCode ierr; 1101 1102 PetscFunctionBegin; 1103 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 1104 if (iascii) { 1105 const TSRosWType rostype; 1106 PetscInt i; 1107 PetscReal abscissa[512]; 1108 char buf[512]; 1109 ierr = TSRosWGetType(ts,&rostype);CHKERRQ(ierr); 1110 ierr = PetscViewerASCIIPrintf(viewer," Rosenbrock-W %s\n",rostype);CHKERRQ(ierr); 1111 ierr = PetscFormatRealArray(buf,sizeof buf,"% 8.6f",tab->s,tab->ASum);CHKERRQ(ierr); 1112 ierr = PetscViewerASCIIPrintf(viewer," Abscissa of A = %s\n",buf);CHKERRQ(ierr); 1113 for (i=0; i<tab->s; i++) abscissa[i] = tab->ASum[i] + tab->Gamma[i]; 1114 ierr = PetscFormatRealArray(buf,sizeof buf,"% 8.6f",tab->s,abscissa);CHKERRQ(ierr); 1115 ierr = PetscViewerASCIIPrintf(viewer," Abscissa of A+Gamma = %s\n",buf);CHKERRQ(ierr); 1116 } 1117 ierr = SNESView(ts->snes,viewer);CHKERRQ(ierr); 1118 PetscFunctionReturn(0); 1119 } 1120 1121 #undef __FUNCT__ 1122 #define __FUNCT__ "TSRosWSetType" 1123 /*@C 1124 TSRosWSetType - Set the type of Rosenbrock-W scheme 1125 1126 Logically collective 1127 1128 Input Parameter: 1129 + ts - timestepping context 1130 - rostype - type of Rosenbrock-W scheme 1131 1132 Level: beginner 1133 1134 .seealso: TSRosWGetType(), TSROSW, TSROSW2M, TSROSW2P, TSROSWRA3PW, TSROSWRA34PW2, TSROSWRODAS3, TSROSWSANDU3, TSROSWASSP3P3S1C, TSROSWLASSP3P4S2C, TSROSWLLSSP3P4S2C, TSROSWARK3 1135 @*/ 1136 PetscErrorCode TSRosWSetType(TS ts,const TSRosWType rostype) 1137 { 1138 PetscErrorCode ierr; 1139 1140 PetscFunctionBegin; 1141 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1142 ierr = PetscTryMethod(ts,"TSRosWSetType_C",(TS,const TSRosWType),(ts,rostype));CHKERRQ(ierr); 1143 PetscFunctionReturn(0); 1144 } 1145 1146 #undef __FUNCT__ 1147 #define __FUNCT__ "TSRosWGetType" 1148 /*@C 1149 TSRosWGetType - Get the type of Rosenbrock-W scheme 1150 1151 Logically collective 1152 1153 Input Parameter: 1154 . ts - timestepping context 1155 1156 Output Parameter: 1157 . rostype - type of Rosenbrock-W scheme 1158 1159 Level: intermediate 1160 1161 .seealso: TSRosWGetType() 1162 @*/ 1163 PetscErrorCode TSRosWGetType(TS ts,const TSRosWType *rostype) 1164 { 1165 PetscErrorCode ierr; 1166 1167 PetscFunctionBegin; 1168 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1169 ierr = PetscUseMethod(ts,"TSRosWGetType_C",(TS,const TSRosWType*),(ts,rostype));CHKERRQ(ierr); 1170 PetscFunctionReturn(0); 1171 } 1172 1173 #undef __FUNCT__ 1174 #define __FUNCT__ "TSRosWSetRecomputeJacobian" 1175 /*@C 1176 TSRosWSetRecomputeJacobian - Set whether to recompute the Jacobian at each stage. The default is to update the Jacobian once per step. 1177 1178 Logically collective 1179 1180 Input Parameter: 1181 + ts - timestepping context 1182 - flg - PETSC_TRUE to recompute the Jacobian at each stage 1183 1184 Level: intermediate 1185 1186 .seealso: TSRosWGetType() 1187 @*/ 1188 PetscErrorCode TSRosWSetRecomputeJacobian(TS ts,PetscBool flg) 1189 { 1190 PetscErrorCode ierr; 1191 1192 PetscFunctionBegin; 1193 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1194 ierr = PetscTryMethod(ts,"TSRosWSetRecomputeJacobian_C",(TS,PetscBool),(ts,flg));CHKERRQ(ierr); 1195 PetscFunctionReturn(0); 1196 } 1197 1198 EXTERN_C_BEGIN 1199 #undef __FUNCT__ 1200 #define __FUNCT__ "TSRosWGetType_RosW" 1201 PetscErrorCode TSRosWGetType_RosW(TS ts,const TSRosWType *rostype) 1202 { 1203 TS_RosW *ros = (TS_RosW*)ts->data; 1204 PetscErrorCode ierr; 1205 1206 PetscFunctionBegin; 1207 if (!ros->tableau) {ierr = TSRosWSetType(ts,TSRosWDefault);CHKERRQ(ierr);} 1208 *rostype = ros->tableau->name; 1209 PetscFunctionReturn(0); 1210 } 1211 #undef __FUNCT__ 1212 #define __FUNCT__ "TSRosWSetType_RosW" 1213 PetscErrorCode TSRosWSetType_RosW(TS ts,const TSRosWType rostype) 1214 { 1215 TS_RosW *ros = (TS_RosW*)ts->data; 1216 PetscErrorCode ierr; 1217 PetscBool match; 1218 RosWTableauLink link; 1219 1220 PetscFunctionBegin; 1221 if (ros->tableau) { 1222 ierr = PetscStrcmp(ros->tableau->name,rostype,&match);CHKERRQ(ierr); 1223 if (match) PetscFunctionReturn(0); 1224 } 1225 for (link = RosWTableauList; link; link=link->next) { 1226 ierr = PetscStrcmp(link->tab.name,rostype,&match);CHKERRQ(ierr); 1227 if (match) { 1228 ierr = TSReset_RosW(ts);CHKERRQ(ierr); 1229 ros->tableau = &link->tab; 1230 PetscFunctionReturn(0); 1231 } 1232 } 1233 SETERRQ1(((PetscObject)ts)->comm,PETSC_ERR_ARG_UNKNOWN_TYPE,"Could not find '%s'",rostype); 1234 PetscFunctionReturn(0); 1235 } 1236 1237 #undef __FUNCT__ 1238 #define __FUNCT__ "TSRosWSetRecomputeJacobian_RosW" 1239 PetscErrorCode TSRosWSetRecomputeJacobian_RosW(TS ts,PetscBool flg) 1240 { 1241 TS_RosW *ros = (TS_RosW*)ts->data; 1242 1243 PetscFunctionBegin; 1244 ros->recompute_jacobian = flg; 1245 PetscFunctionReturn(0); 1246 } 1247 EXTERN_C_END 1248 1249 /* ------------------------------------------------------------ */ 1250 /*MC 1251 TSROSW - ODE solver using Rosenbrock-W schemes 1252 1253 These methods are intended for problems with well-separated time scales, especially when a slow scale is strongly 1254 nonlinear such that it is expensive to solve with a fully implicit method. The user should provide the stiff part 1255 of the equation using TSSetIFunction() and the non-stiff part with TSSetRHSFunction(). 1256 1257 Notes: 1258 This method currently only works with autonomous ODE and DAE. 1259 1260 Developer notes: 1261 Rosenbrock-W methods are typically specified for autonomous ODE 1262 1263 $ xdot = f(x) 1264 1265 by the stage equations 1266 1267 $ k_i = h f(x_0 + sum_j alpha_ij k_j) + h J sum_j gamma_ij k_j 1268 1269 and step completion formula 1270 1271 $ x_1 = x_0 + sum_j b_j k_j 1272 1273 with step size h and coefficients alpha_ij, gamma_ij, and b_i. Implementing the method in this form would require f(x) 1274 and the Jacobian J to be available, in addition to the shifted matrix I - h gamma_ii J. Following Hairer and Wanner, 1275 we define new variables for the stage equations 1276 1277 $ y_i = gamma_ij k_j 1278 1279 The k_j can be recovered because Gamma is invertible. Let C be the lower triangular part of Gamma^{-1} and define 1280 1281 $ A = Alpha Gamma^{-1}, bt^T = b^T Gamma^{-i} 1282 1283 to rewrite the method as 1284 1285 $ [M/(h gamma_ii) - J] y_i = f(x_0 + sum_j a_ij y_j) + M sum_j (c_ij/h) y_j 1286 $ x_1 = x_0 + sum_j bt_j y_j 1287 1288 where we have introduced the mass matrix M. Continue by defining 1289 1290 $ ydot_i = 1/(h gamma_ii) y_i - sum_j (c_ij/h) y_j 1291 1292 or, more compactly in tensor notation 1293 1294 $ Ydot = 1/h (Gamma^{-1} \otimes I) Y . 1295 1296 Note that Gamma^{-1} is lower triangular. With this definition of Ydot in terms of known quantities and the current 1297 stage y_i, the stage equations reduce to performing one Newton step (typically with a lagged Jacobian) on the 1298 equation 1299 1300 $ g(x_0 + sum_j a_ij y_j + y_i, ydot_i) = 0 1301 1302 with initial guess y_i = 0. 1303 1304 Level: beginner 1305 1306 .seealso: TSCreate(), TS, TSSetType(), TSRosWSetType(), TSRosWRegister() 1307 1308 M*/ 1309 EXTERN_C_BEGIN 1310 #undef __FUNCT__ 1311 #define __FUNCT__ "TSCreate_RosW" 1312 PetscErrorCode TSCreate_RosW(TS ts) 1313 { 1314 TS_RosW *ros; 1315 PetscErrorCode ierr; 1316 1317 PetscFunctionBegin; 1318 #if !defined(PETSC_USE_DYNAMIC_LIBRARIES) 1319 ierr = TSRosWInitializePackage(PETSC_NULL);CHKERRQ(ierr); 1320 #endif 1321 1322 ts->ops->reset = TSReset_RosW; 1323 ts->ops->destroy = TSDestroy_RosW; 1324 ts->ops->view = TSView_RosW; 1325 ts->ops->setup = TSSetUp_RosW; 1326 ts->ops->step = TSStep_RosW; 1327 ts->ops->interpolate = TSInterpolate_RosW; 1328 ts->ops->evaluatestep = TSEvaluateStep_RosW; 1329 ts->ops->setfromoptions = TSSetFromOptions_RosW; 1330 ts->ops->snesfunction = SNESTSFormFunction_RosW; 1331 ts->ops->snesjacobian = SNESTSFormJacobian_RosW; 1332 1333 ierr = PetscNewLog(ts,TS_RosW,&ros);CHKERRQ(ierr); 1334 ts->data = (void*)ros; 1335 1336 ierr = PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSRosWGetType_C","TSRosWGetType_RosW",TSRosWGetType_RosW);CHKERRQ(ierr); 1337 ierr = PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSRosWSetType_C","TSRosWSetType_RosW",TSRosWSetType_RosW);CHKERRQ(ierr); 1338 ierr = PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSRosWSetRecomputeJacobian_C","TSRosWSetRecomputeJacobian_RosW",TSRosWSetRecomputeJacobian_RosW);CHKERRQ(ierr); 1339 PetscFunctionReturn(0); 1340 } 1341 EXTERN_C_END 1342