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 = VecCopy(ts->vec_sol,ros->VecSolPrev);CHKERRQ(ierr);/*move this at the end*/ 801 for (i=0; i<s; i++) { 802 ros->stage_time = ts->ptime + h*ASum[i]; 803 if (GammaZeroDiag[i]) { 804 ros->stage_explicit = PETSC_TRUE; 805 ros->shift = 1./h; 806 } else { 807 ros->stage_explicit = PETSC_FALSE; 808 ros->shift = 1./(h*Gamma[i*s+i]); 809 } 810 811 ierr = VecCopy(ts->vec_sol,Zstage);CHKERRQ(ierr); 812 for (j=0; j<i; j++) w[j] = At[i*s+j]; 813 ierr = VecMAXPY(Zstage,i,w,Y);CHKERRQ(ierr); 814 815 for (j=0; j<i; j++) w[j] = 1./h * GammaInv[i*s+j]; 816 ierr = VecZeroEntries(Zdot);CHKERRQ(ierr); 817 ierr = VecMAXPY(Zdot,i,w,Y);CHKERRQ(ierr); 818 819 /* Initial guess taken from last stage */ 820 ierr = VecZeroEntries(Y[i]);CHKERRQ(ierr); 821 822 if (!ros->stage_explicit) { 823 if (!ros->recompute_jacobian && !i) { 824 ierr = SNESSetLagJacobian(snes,-2);CHKERRQ(ierr); /* Recompute the Jacobian on this solve, but not again */ 825 } 826 ierr = SNESSolve(snes,PETSC_NULL,Y[i]);CHKERRQ(ierr); 827 ierr = SNESGetIterationNumber(snes,&its);CHKERRQ(ierr); 828 ierr = SNESGetLinearSolveIterations(snes,&lits);CHKERRQ(ierr); 829 ts->snes_its += its; ts->ksp_its += lits; 830 ierr = TSGetAdapt(ts,&adapt);CHKERRQ(ierr); 831 ierr = TSAdaptCheckStage(adapt,ts,&accept);CHKERRQ(ierr); 832 if (!accept) goto reject_step; 833 } else { 834 Mat J,Jp; 835 ierr = VecZeroEntries(Ydot);CHKERRQ(ierr); /* Evaluate Y[i]=G(t,Ydot=0,Zstage) */ 836 ierr = TSComputeIFunction(ts,ros->stage_time,Zstage,Ydot,Y[i],PETSC_FALSE);CHKERRQ(ierr); 837 ierr = VecScale(Y[i],-1.0); 838 ierr = VecAXPY(Y[i],-1.0,Zdot);CHKERRQ(ierr); /*Y[i]=F(Zstage)-Zdot[=GammaInv*Y]*/ 839 840 ierr = VecZeroEntries(Zstage);CHKERRQ(ierr); /* Zstage = GammaExplicitCorr[i,j] * Y[j] */ 841 for (j=0; j<i; j++) w[j] = GammaExplicitCorr[i*s+j]; 842 ierr = VecMAXPY(Zstage,i,w,Y);CHKERRQ(ierr); 843 /*Y[i] += Y[i] + Jac*Zstage[=Jac*GammaExplicitCorr[i,j] * Y[j]] */ 844 str = SAME_NONZERO_PATTERN; 845 ierr = TSGetIJacobian(ts,&J,&Jp,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr); 846 ierr = TSComputeIJacobian(ts,ros->stage_time,ts->vec_sol,Ydot,0,&J,&Jp,&str,PETSC_FALSE);CHKERRQ(ierr); 847 ierr = MatMult(J,Zstage,Zdot); 848 849 ierr = VecAXPY(Y[i],-1.0,Zdot);CHKERRQ(ierr); 850 ierr = VecScale(Y[i],h); 851 ts->ksp_its += 1; 852 } 853 } 854 ierr = TSEvaluateStep(ts,tab->order,ts->vec_sol,PETSC_NULL);CHKERRQ(ierr); 855 ros->status = TS_STEP_PENDING; 856 857 /* Register only the current method as a candidate because we're not supporting multiple candidates yet. */ 858 ierr = TSGetAdapt(ts,&adapt);CHKERRQ(ierr); 859 ierr = TSAdaptCandidatesClear(adapt);CHKERRQ(ierr); 860 ierr = TSAdaptCandidateAdd(adapt,tab->name,tab->order,1,tab->ccfl,1.*tab->s,PETSC_TRUE);CHKERRQ(ierr); 861 ierr = TSAdaptChoose(adapt,ts,ts->time_step,&next_scheme,&next_time_step,&accept);CHKERRQ(ierr); 862 if (accept) { 863 /* ignore next_scheme for now */ 864 ts->ptime += ts->time_step; 865 ts->time_step = next_time_step; 866 ts->steps++; 867 ros->status = TS_STEP_COMPLETE; 868 break; 869 } else { /* Roll back the current step */ 870 for (i=0; i<s; i++) w[i] = -tab->bt[i]; 871 ierr = VecMAXPY(ts->vec_sol,s,w,Y);CHKERRQ(ierr); 872 ts->time_step = next_time_step; 873 ros->status = TS_STEP_INCOMPLETE; 874 } 875 reject_step: continue; 876 } 877 if (ros->status != TS_STEP_COMPLETE && !ts->reason) ts->reason = TS_DIVERGED_STEP_REJECTED; 878 PetscFunctionReturn(0); 879 } 880 881 #undef __FUNCT__ 882 #define __FUNCT__ "TSInterpolate_RosW" 883 static PetscErrorCode TSInterpolate_RosW(TS ts,PetscReal itime,Vec X) 884 { 885 TS_RosW *ros = (TS_RosW*)ts->data; 886 PetscInt s = ros->tableau->s,pinterp = ros->tableau->pinterp,i,j; 887 PetscReal h; 888 PetscReal tt,t; 889 PetscScalar *bt; 890 const PetscReal *Bt = ros->tableau->binterpt; 891 PetscErrorCode ierr; 892 const PetscReal *GammaInv = ros->tableau->GammaInv; 893 PetscScalar *w = ros->work; 894 Vec *Y = ros->Y; 895 896 PetscFunctionBegin; 897 if (!Bt) SETERRQ1(((PetscObject)ts)->comm,PETSC_ERR_SUP,"TSRosW %s does not have an interpolation formula",ros->tableau->name); 898 899 switch (ros->status) { 900 case TS_STEP_INCOMPLETE: 901 case TS_STEP_PENDING: 902 h = ts->time_step; 903 t = (itime - ts->ptime)/h; 904 break; 905 case TS_STEP_COMPLETE: 906 h = ts->time_step_prev; 907 t = (itime - ts->ptime)/h + 1; /* In the interval [0,1] */ 908 break; 909 default: SETERRQ(((PetscObject)ts)->comm,PETSC_ERR_PLIB,"Invalid TSStepStatus"); 910 } 911 ierr = PetscMalloc(s*sizeof(bt[0]),&bt);CHKERRQ(ierr); 912 for (i=0; i<s; i++) bt[i] = 0; 913 for (j=0,tt=t; j<pinterp; j++,tt*=t) { 914 for (i=0; i<s; i++) { 915 bt[i] += Bt[i*pinterp+j] * tt; 916 } 917 } 918 919 /* y(t+tt*h) = y(t) + Sum bt(tt) * GammaInv * Ydot */ 920 /*X<-0*/ 921 ierr = VecZeroEntries(X);CHKERRQ(ierr); 922 923 /*X<- Sum bt_i * GammaInv(i,1:i) * Y(1:i) */ 924 for (j=0; j<s; j++) w[j]=0; 925 for (j=0; j<s; j++) { 926 for (i=j; i<s; i++) { 927 w[j] += bt[i]*GammaInv[i*s+j]; 928 } 929 } 930 ierr = VecMAXPY(X,i,w,Y);CHKERRQ(ierr); 931 932 /*X<-y(t) + X*/ 933 ierr = VecAXPY(X,1.0,ros->VecSolPrev);CHKERRQ(ierr); 934 935 ierr = PetscFree(bt);CHKERRQ(ierr); 936 937 PetscFunctionReturn(0); 938 } 939 940 /*------------------------------------------------------------*/ 941 #undef __FUNCT__ 942 #define __FUNCT__ "TSReset_RosW" 943 static PetscErrorCode TSReset_RosW(TS ts) 944 { 945 TS_RosW *ros = (TS_RosW*)ts->data; 946 PetscInt s; 947 PetscErrorCode ierr; 948 949 PetscFunctionBegin; 950 if (!ros->tableau) PetscFunctionReturn(0); 951 s = ros->tableau->s; 952 ierr = VecDestroyVecs(s,&ros->Y);CHKERRQ(ierr); 953 ierr = VecDestroy(&ros->Ydot);CHKERRQ(ierr); 954 ierr = VecDestroy(&ros->Ystage);CHKERRQ(ierr); 955 ierr = VecDestroy(&ros->Zdot);CHKERRQ(ierr); 956 ierr = VecDestroy(&ros->Zstage);CHKERRQ(ierr); 957 ierr = VecDestroy(&ros->VecSolPrev);CHKERRQ(ierr); 958 ierr = PetscFree(ros->work);CHKERRQ(ierr); 959 PetscFunctionReturn(0); 960 } 961 962 #undef __FUNCT__ 963 #define __FUNCT__ "TSDestroy_RosW" 964 static PetscErrorCode TSDestroy_RosW(TS ts) 965 { 966 PetscErrorCode ierr; 967 968 PetscFunctionBegin; 969 ierr = TSReset_RosW(ts);CHKERRQ(ierr); 970 ierr = PetscFree(ts->data);CHKERRQ(ierr); 971 ierr = PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSRosWGetType_C","",PETSC_NULL);CHKERRQ(ierr); 972 ierr = PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSRosWSetType_C","",PETSC_NULL);CHKERRQ(ierr); 973 ierr = PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSRosWSetRecomputeJacobian_C","",PETSC_NULL);CHKERRQ(ierr); 974 PetscFunctionReturn(0); 975 } 976 977 /* 978 This defines the nonlinear equation that is to be solved with SNES 979 G(U) = F[t0+Theta*dt, U, (U-U0)*shift] = 0 980 */ 981 #undef __FUNCT__ 982 #define __FUNCT__ "SNESTSFormFunction_RosW" 983 static PetscErrorCode SNESTSFormFunction_RosW(SNES snes,Vec X,Vec F,TS ts) 984 { 985 TS_RosW *ros = (TS_RosW*)ts->data; 986 PetscErrorCode ierr; 987 988 PetscFunctionBegin; 989 ierr = VecWAXPY(ros->Ydot,ros->shift,X,ros->Zdot);CHKERRQ(ierr); /* Ydot = shift*X + Zdot */ 990 ierr = VecWAXPY(ros->Ystage,1.0,X,ros->Zstage);CHKERRQ(ierr); /* Ystage = X + Zstage */ 991 ierr = TSComputeIFunction(ts,ros->stage_time,ros->Ystage,ros->Ydot,F,PETSC_FALSE);CHKERRQ(ierr); 992 PetscFunctionReturn(0); 993 } 994 995 #undef __FUNCT__ 996 #define __FUNCT__ "SNESTSFormJacobian_RosW" 997 static PetscErrorCode SNESTSFormJacobian_RosW(SNES snes,Vec X,Mat *A,Mat *B,MatStructure *str,TS ts) 998 { 999 TS_RosW *ros = (TS_RosW*)ts->data; 1000 PetscErrorCode ierr; 1001 1002 PetscFunctionBegin; 1003 /* ros->Ydot and ros->Ystage have already been computed in SNESTSFormFunction_RosW (SNES guarantees this) */ 1004 ierr = TSComputeIJacobian(ts,ros->stage_time,ros->Ystage,ros->Ydot,ros->shift,A,B,str,PETSC_TRUE);CHKERRQ(ierr); 1005 PetscFunctionReturn(0); 1006 } 1007 1008 #undef __FUNCT__ 1009 #define __FUNCT__ "TSSetUp_RosW" 1010 static PetscErrorCode TSSetUp_RosW(TS ts) 1011 { 1012 TS_RosW *ros = (TS_RosW*)ts->data; 1013 RosWTableau tab = ros->tableau; 1014 PetscInt s = tab->s; 1015 PetscErrorCode ierr; 1016 1017 PetscFunctionBegin; 1018 if (!ros->tableau) { 1019 ierr = TSRosWSetType(ts,TSRosWDefault);CHKERRQ(ierr); 1020 } 1021 ierr = VecDuplicateVecs(ts->vec_sol,s,&ros->Y);CHKERRQ(ierr); 1022 ierr = VecDuplicate(ts->vec_sol,&ros->Ydot);CHKERRQ(ierr); 1023 ierr = VecDuplicate(ts->vec_sol,&ros->Ystage);CHKERRQ(ierr); 1024 ierr = VecDuplicate(ts->vec_sol,&ros->Zdot);CHKERRQ(ierr); 1025 ierr = VecDuplicate(ts->vec_sol,&ros->Zstage);CHKERRQ(ierr); 1026 ierr = VecDuplicate(ts->vec_sol,&ros->VecSolPrev);CHKERRQ(ierr); 1027 ierr = PetscMalloc(s*sizeof(ros->work[0]),&ros->work);CHKERRQ(ierr); 1028 PetscFunctionReturn(0); 1029 } 1030 /*------------------------------------------------------------*/ 1031 1032 #undef __FUNCT__ 1033 #define __FUNCT__ "TSSetFromOptions_RosW" 1034 static PetscErrorCode TSSetFromOptions_RosW(TS ts) 1035 { 1036 TS_RosW *ros = (TS_RosW*)ts->data; 1037 PetscErrorCode ierr; 1038 char rostype[256]; 1039 1040 PetscFunctionBegin; 1041 ierr = PetscOptionsHead("RosW ODE solver options");CHKERRQ(ierr); 1042 { 1043 RosWTableauLink link; 1044 PetscInt count,choice; 1045 PetscBool flg; 1046 const char **namelist; 1047 SNES snes; 1048 1049 ierr = PetscStrncpy(rostype,TSRosWDefault,sizeof rostype);CHKERRQ(ierr); 1050 for (link=RosWTableauList,count=0; link; link=link->next,count++) ; 1051 ierr = PetscMalloc(count*sizeof(char*),&namelist);CHKERRQ(ierr); 1052 for (link=RosWTableauList,count=0; link; link=link->next,count++) namelist[count] = link->tab.name; 1053 ierr = PetscOptionsEList("-ts_rosw_type","Family of Rosenbrock-W method","TSRosWSetType",(const char*const*)namelist,count,rostype,&choice,&flg);CHKERRQ(ierr); 1054 ierr = TSRosWSetType(ts,flg ? namelist[choice] : rostype);CHKERRQ(ierr); 1055 ierr = PetscFree(namelist);CHKERRQ(ierr); 1056 1057 ierr = PetscOptionsBool("-ts_rosw_recompute_jacobian","Recompute the Jacobian at each stage","TSRosWSetRecomputeJacobian",ros->recompute_jacobian,&ros->recompute_jacobian,PETSC_NULL);CHKERRQ(ierr); 1058 1059 /* Rosenbrock methods are linearly implicit, so set that unless the user has specifically asked for something else */ 1060 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 1061 if (!((PetscObject)snes)->type_name) { 1062 ierr = SNESSetType(snes,SNESKSPONLY);CHKERRQ(ierr); 1063 } 1064 ierr = SNESSetFromOptions(snes);CHKERRQ(ierr); 1065 } 1066 ierr = PetscOptionsTail();CHKERRQ(ierr); 1067 PetscFunctionReturn(0); 1068 } 1069 1070 #undef __FUNCT__ 1071 #define __FUNCT__ "PetscFormatRealArray" 1072 static PetscErrorCode PetscFormatRealArray(char buf[],size_t len,const char *fmt,PetscInt n,const PetscReal x[]) 1073 { 1074 PetscErrorCode ierr; 1075 PetscInt i; 1076 size_t left,count; 1077 char *p; 1078 1079 PetscFunctionBegin; 1080 for (i=0,p=buf,left=len; i<n; i++) { 1081 ierr = PetscSNPrintfCount(p,left,fmt,&count,x[i]);CHKERRQ(ierr); 1082 if (count >= left) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Insufficient space in buffer"); 1083 left -= count; 1084 p += count; 1085 *p++ = ' '; 1086 } 1087 p[i ? 0 : -1] = 0; 1088 PetscFunctionReturn(0); 1089 } 1090 1091 #undef __FUNCT__ 1092 #define __FUNCT__ "TSView_RosW" 1093 static PetscErrorCode TSView_RosW(TS ts,PetscViewer viewer) 1094 { 1095 TS_RosW *ros = (TS_RosW*)ts->data; 1096 RosWTableau tab = ros->tableau; 1097 PetscBool iascii; 1098 PetscErrorCode ierr; 1099 1100 PetscFunctionBegin; 1101 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 1102 if (iascii) { 1103 const TSRosWType rostype; 1104 PetscInt i; 1105 PetscReal abscissa[512]; 1106 char buf[512]; 1107 ierr = TSRosWGetType(ts,&rostype);CHKERRQ(ierr); 1108 ierr = PetscViewerASCIIPrintf(viewer," Rosenbrock-W %s\n",rostype);CHKERRQ(ierr); 1109 ierr = PetscFormatRealArray(buf,sizeof buf,"% 8.6f",tab->s,tab->ASum);CHKERRQ(ierr); 1110 ierr = PetscViewerASCIIPrintf(viewer," Abscissa of A = %s\n",buf);CHKERRQ(ierr); 1111 for (i=0; i<tab->s; i++) abscissa[i] = tab->ASum[i] + tab->Gamma[i]; 1112 ierr = PetscFormatRealArray(buf,sizeof buf,"% 8.6f",tab->s,abscissa);CHKERRQ(ierr); 1113 ierr = PetscViewerASCIIPrintf(viewer," Abscissa of A+Gamma = %s\n",buf);CHKERRQ(ierr); 1114 } 1115 ierr = SNESView(ts->snes,viewer);CHKERRQ(ierr); 1116 PetscFunctionReturn(0); 1117 } 1118 1119 #undef __FUNCT__ 1120 #define __FUNCT__ "TSRosWSetType" 1121 /*@C 1122 TSRosWSetType - Set the type of Rosenbrock-W scheme 1123 1124 Logically collective 1125 1126 Input Parameter: 1127 + ts - timestepping context 1128 - rostype - type of Rosenbrock-W scheme 1129 1130 Level: beginner 1131 1132 .seealso: TSRosWGetType(), TSROSW, TSROSW2M, TSROSW2P, TSROSWRA3PW, TSROSWRA34PW2, TSROSWRODAS3, TSROSWSANDU3, TSROSWASSP3P3S1C, TSROSWLASSP3P4S2C, TSROSWLLSSP3P4S2C, TSROSWARK3 1133 @*/ 1134 PetscErrorCode TSRosWSetType(TS ts,const TSRosWType rostype) 1135 { 1136 PetscErrorCode ierr; 1137 1138 PetscFunctionBegin; 1139 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1140 ierr = PetscTryMethod(ts,"TSRosWSetType_C",(TS,const TSRosWType),(ts,rostype));CHKERRQ(ierr); 1141 PetscFunctionReturn(0); 1142 } 1143 1144 #undef __FUNCT__ 1145 #define __FUNCT__ "TSRosWGetType" 1146 /*@C 1147 TSRosWGetType - Get the type of Rosenbrock-W scheme 1148 1149 Logically collective 1150 1151 Input Parameter: 1152 . ts - timestepping context 1153 1154 Output Parameter: 1155 . rostype - type of Rosenbrock-W scheme 1156 1157 Level: intermediate 1158 1159 .seealso: TSRosWGetType() 1160 @*/ 1161 PetscErrorCode TSRosWGetType(TS ts,const TSRosWType *rostype) 1162 { 1163 PetscErrorCode ierr; 1164 1165 PetscFunctionBegin; 1166 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1167 ierr = PetscUseMethod(ts,"TSRosWGetType_C",(TS,const TSRosWType*),(ts,rostype));CHKERRQ(ierr); 1168 PetscFunctionReturn(0); 1169 } 1170 1171 #undef __FUNCT__ 1172 #define __FUNCT__ "TSRosWSetRecomputeJacobian" 1173 /*@C 1174 TSRosWSetRecomputeJacobian - Set whether to recompute the Jacobian at each stage. The default is to update the Jacobian once per step. 1175 1176 Logically collective 1177 1178 Input Parameter: 1179 + ts - timestepping context 1180 - flg - PETSC_TRUE to recompute the Jacobian at each stage 1181 1182 Level: intermediate 1183 1184 .seealso: TSRosWGetType() 1185 @*/ 1186 PetscErrorCode TSRosWSetRecomputeJacobian(TS ts,PetscBool flg) 1187 { 1188 PetscErrorCode ierr; 1189 1190 PetscFunctionBegin; 1191 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1192 ierr = PetscTryMethod(ts,"TSRosWSetRecomputeJacobian_C",(TS,PetscBool),(ts,flg));CHKERRQ(ierr); 1193 PetscFunctionReturn(0); 1194 } 1195 1196 EXTERN_C_BEGIN 1197 #undef __FUNCT__ 1198 #define __FUNCT__ "TSRosWGetType_RosW" 1199 PetscErrorCode TSRosWGetType_RosW(TS ts,const TSRosWType *rostype) 1200 { 1201 TS_RosW *ros = (TS_RosW*)ts->data; 1202 PetscErrorCode ierr; 1203 1204 PetscFunctionBegin; 1205 if (!ros->tableau) {ierr = TSRosWSetType(ts,TSRosWDefault);CHKERRQ(ierr);} 1206 *rostype = ros->tableau->name; 1207 PetscFunctionReturn(0); 1208 } 1209 #undef __FUNCT__ 1210 #define __FUNCT__ "TSRosWSetType_RosW" 1211 PetscErrorCode TSRosWSetType_RosW(TS ts,const TSRosWType rostype) 1212 { 1213 TS_RosW *ros = (TS_RosW*)ts->data; 1214 PetscErrorCode ierr; 1215 PetscBool match; 1216 RosWTableauLink link; 1217 1218 PetscFunctionBegin; 1219 if (ros->tableau) { 1220 ierr = PetscStrcmp(ros->tableau->name,rostype,&match);CHKERRQ(ierr); 1221 if (match) PetscFunctionReturn(0); 1222 } 1223 for (link = RosWTableauList; link; link=link->next) { 1224 ierr = PetscStrcmp(link->tab.name,rostype,&match);CHKERRQ(ierr); 1225 if (match) { 1226 ierr = TSReset_RosW(ts);CHKERRQ(ierr); 1227 ros->tableau = &link->tab; 1228 PetscFunctionReturn(0); 1229 } 1230 } 1231 SETERRQ1(((PetscObject)ts)->comm,PETSC_ERR_ARG_UNKNOWN_TYPE,"Could not find '%s'",rostype); 1232 PetscFunctionReturn(0); 1233 } 1234 1235 #undef __FUNCT__ 1236 #define __FUNCT__ "TSRosWSetRecomputeJacobian_RosW" 1237 PetscErrorCode TSRosWSetRecomputeJacobian_RosW(TS ts,PetscBool flg) 1238 { 1239 TS_RosW *ros = (TS_RosW*)ts->data; 1240 1241 PetscFunctionBegin; 1242 ros->recompute_jacobian = flg; 1243 PetscFunctionReturn(0); 1244 } 1245 EXTERN_C_END 1246 1247 /* ------------------------------------------------------------ */ 1248 /*MC 1249 TSROSW - ODE solver using Rosenbrock-W schemes 1250 1251 These methods are intended for problems with well-separated time scales, especially when a slow scale is strongly 1252 nonlinear such that it is expensive to solve with a fully implicit method. The user should provide the stiff part 1253 of the equation using TSSetIFunction() and the non-stiff part with TSSetRHSFunction(). 1254 1255 Notes: 1256 This method currently only works with autonomous ODE and DAE. 1257 1258 Developer notes: 1259 Rosenbrock-W methods are typically specified for autonomous ODE 1260 1261 $ xdot = f(x) 1262 1263 by the stage equations 1264 1265 $ k_i = h f(x_0 + sum_j alpha_ij k_j) + h J sum_j gamma_ij k_j 1266 1267 and step completion formula 1268 1269 $ x_1 = x_0 + sum_j b_j k_j 1270 1271 with step size h and coefficients alpha_ij, gamma_ij, and b_i. Implementing the method in this form would require f(x) 1272 and the Jacobian J to be available, in addition to the shifted matrix I - h gamma_ii J. Following Hairer and Wanner, 1273 we define new variables for the stage equations 1274 1275 $ y_i = gamma_ij k_j 1276 1277 The k_j can be recovered because Gamma is invertible. Let C be the lower triangular part of Gamma^{-1} and define 1278 1279 $ A = Alpha Gamma^{-1}, bt^T = b^T Gamma^{-i} 1280 1281 to rewrite the method as 1282 1283 $ [M/(h gamma_ii) - J] y_i = f(x_0 + sum_j a_ij y_j) + M sum_j (c_ij/h) y_j 1284 $ x_1 = x_0 + sum_j bt_j y_j 1285 1286 where we have introduced the mass matrix M. Continue by defining 1287 1288 $ ydot_i = 1/(h gamma_ii) y_i - sum_j (c_ij/h) y_j 1289 1290 or, more compactly in tensor notation 1291 1292 $ Ydot = 1/h (Gamma^{-1} \otimes I) Y . 1293 1294 Note that Gamma^{-1} is lower triangular. With this definition of Ydot in terms of known quantities and the current 1295 stage y_i, the stage equations reduce to performing one Newton step (typically with a lagged Jacobian) on the 1296 equation 1297 1298 $ g(x_0 + sum_j a_ij y_j + y_i, ydot_i) = 0 1299 1300 with initial guess y_i = 0. 1301 1302 Level: beginner 1303 1304 .seealso: TSCreate(), TS, TSSetType(), TSRosWSetType(), TSRosWRegister() 1305 1306 M*/ 1307 EXTERN_C_BEGIN 1308 #undef __FUNCT__ 1309 #define __FUNCT__ "TSCreate_RosW" 1310 PetscErrorCode TSCreate_RosW(TS ts) 1311 { 1312 TS_RosW *ros; 1313 PetscErrorCode ierr; 1314 1315 PetscFunctionBegin; 1316 #if !defined(PETSC_USE_DYNAMIC_LIBRARIES) 1317 ierr = TSRosWInitializePackage(PETSC_NULL);CHKERRQ(ierr); 1318 #endif 1319 1320 ts->ops->reset = TSReset_RosW; 1321 ts->ops->destroy = TSDestroy_RosW; 1322 ts->ops->view = TSView_RosW; 1323 ts->ops->setup = TSSetUp_RosW; 1324 ts->ops->step = TSStep_RosW; 1325 ts->ops->interpolate = TSInterpolate_RosW; 1326 ts->ops->evaluatestep = TSEvaluateStep_RosW; 1327 ts->ops->setfromoptions = TSSetFromOptions_RosW; 1328 ts->ops->snesfunction = SNESTSFormFunction_RosW; 1329 ts->ops->snesjacobian = SNESTSFormJacobian_RosW; 1330 1331 ierr = PetscNewLog(ts,TS_RosW,&ros);CHKERRQ(ierr); 1332 ts->data = (void*)ros; 1333 1334 ierr = PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSRosWGetType_C","TSRosWGetType_RosW",TSRosWGetType_RosW);CHKERRQ(ierr); 1335 ierr = PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSRosWSetType_C","TSRosWSetType_RosW",TSRosWSetType_RosW);CHKERRQ(ierr); 1336 ierr = PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSRosWSetRecomputeJacobian_C","TSRosWSetRecomputeJacobian_RosW",TSRosWSetRecomputeJacobian_RosW);CHKERRQ(ierr); 1337 PetscFunctionReturn(0); 1338 } 1339 EXTERN_C_END 1340