1 /* 2 Code for time stepping with the General Linear with Error Estimation method 3 4 5 Notes: 6 The general system is written as 7 8 Udot = F(t,U) 9 10 */ 11 #include <petsc/private/tsimpl.h> /*I "petscts.h" I*/ 12 #include <petscdm.h> 13 14 static PetscBool cited = PETSC_FALSE; 15 static const char citation[] = 16 "@ARTICLE{Constantinescu_TR2016b,\n" 17 " author = {Constantinescu, E.M.},\n" 18 " title = {Estimating Global Errors in Time Stepping},\n" 19 " journal = {ArXiv e-prints},\n" 20 " year = 2016,\n" 21 " adsurl = {http://adsabs.harvard.edu/abs/2015arXiv150305166C}\n}\n"; 22 23 static TSGLEEType TSGLEEDefaultType = TSGLEE35; 24 static PetscBool TSGLEERegisterAllCalled; 25 static PetscBool TSGLEEPackageInitialized; 26 static PetscInt explicit_stage_time_id; 27 28 typedef struct _GLEETableau *GLEETableau; 29 struct _GLEETableau { 30 char *name; 31 PetscInt order; /* Classical approximation order of the method i*/ 32 PetscInt s; /* Number of stages */ 33 PetscInt r; /* Number of steps */ 34 PetscReal gamma; /* LTE ratio */ 35 PetscReal *A,*B,*U,*V,*S,*F,*c; /* Tableau */ 36 PetscReal *Fembed; /* Embedded final method coefficients */ 37 PetscReal *Ferror; /* Coefficients for computing error */ 38 PetscReal *Serror; /* Coefficients for initializing the error */ 39 PetscInt pinterp; /* Interpolation order */ 40 PetscReal *binterp; /* Interpolation coefficients */ 41 PetscReal ccfl; /* Placeholder for CFL coefficient relative to forward Euler */ 42 }; 43 typedef struct _GLEETableauLink *GLEETableauLink; 44 struct _GLEETableauLink { 45 struct _GLEETableau tab; 46 GLEETableauLink next; 47 }; 48 static GLEETableauLink GLEETableauList; 49 50 typedef struct { 51 GLEETableau tableau; 52 Vec *Y; /* Solution vector (along with auxiliary solution y~ or eps) */ 53 Vec *X; /* Temporary solution vector */ 54 Vec *YStage; /* Stage values */ 55 Vec *YdotStage; /* Stage right hand side */ 56 Vec W; /* Right-hand-side for implicit stage solve */ 57 Vec Ydot; /* Work vector holding Ydot during residual evaluation */ 58 Vec yGErr; /* Vector holding the global error after a step is completed */ 59 PetscScalar *work; /* Scalar work */ 60 PetscReal scoeff; /* shift = scoeff/dt */ 61 PetscReal stage_time; 62 TSStepStatus status; 63 } TS_GLEE; 64 65 /*MC 66 TSGLEE23 - Second order three stage GLEE method 67 68 This method has three stages. 69 s = 3, r = 2 70 71 Level: advanced 72 73 .seealso: TSGLEE 74 */ 75 /*MC 76 TSGLEE24 - Second order four stage GLEE method 77 78 This method has four stages. 79 s = 4, r = 2 80 81 Level: advanced 82 83 .seealso: TSGLEE 84 */ 85 /*MC 86 TSGLEE25i - Second order five stage GLEE method 87 88 This method has five stages. 89 s = 5, r = 2 90 91 Level: advanced 92 93 .seealso: TSGLEE 94 */ 95 /*MC 96 TSGLEE35 - Third order five stage GLEE method 97 98 This method has five stages. 99 s = 5, r = 2 100 101 Level: advanced 102 103 .seealso: TSGLEE 104 */ 105 /*MC 106 TSGLEEEXRK2A - Second order six stage GLEE method 107 108 This method has six stages. 109 s = 6, r = 2 110 111 Level: advanced 112 113 .seealso: TSGLEE 114 */ 115 /*MC 116 TSGLEERK32G1 - Third order eight stage GLEE method 117 118 This method has eight stages. 119 s = 8, r = 2 120 121 Level: advanced 122 123 .seealso: TSGLEE 124 */ 125 /*MC 126 TSGLEERK285EX - Second order nine stage GLEE method 127 128 This method has nine stages. 129 s = 9, r = 2 130 131 Level: advanced 132 133 .seealso: TSGLEE 134 */ 135 136 #undef __FUNCT__ 137 #define __FUNCT__ "TSGLEERegisterAll" 138 /*@C 139 TSGLEERegisterAll - Registers all of the General Linear with Error Estimation methods in TSGLEE 140 141 Not Collective, but should be called by all processes which will need the schemes to be registered 142 143 Level: advanced 144 145 .keywords: TS, TSGLEE, register, all 146 147 .seealso: TSGLEERegisterDestroy() 148 @*/ 149 PetscErrorCode TSGLEERegisterAll(void) 150 { 151 PetscErrorCode ierr; 152 153 PetscFunctionBegin; 154 if (TSGLEERegisterAllCalled) PetscFunctionReturn(0); 155 TSGLEERegisterAllCalled = PETSC_TRUE; 156 157 { 158 /* y-eps form */ 159 const PetscInt 160 p = 1, 161 s = 3, 162 r = 2; 163 const PetscReal 164 gamma = 0.5, 165 A[3][3] = {{1.0,0,0},{0,0.5,0},{0,0.5,0.5}}, 166 B[2][3] = {{1.0,0,0},{-2.0,1.0,1.0}}, 167 U[3][2] = {{1.0,0},{1.0,0.5},{1.0,0.5}}, 168 V[2][2] = {{1,0},{0,1}}, 169 S[2] = {1,0}, 170 F[2] = {1,0}, 171 Fembed[2] = {1,1-gamma}, 172 Ferror[2] = {0,1}, 173 Serror[2] = {1,0}; 174 ierr = TSGLEERegister(TSGLEEi1,p,s,r,gamma,&A[0][0],&B[0][0],&U[0][0],&V[0][0],S,F,NULL,Fembed,Ferror,Serror,0,NULL);CHKERRQ(ierr); 175 } 176 { 177 /* y-eps form */ 178 const PetscInt 179 p = 2, 180 s = 3, 181 r = 2; 182 const PetscReal 183 gamma = 0, 184 A[3][3] = {{0,0,0},{1,0,0},{0.25,0.25,0}}, 185 B[2][3] = {{1.0/12.0,1.0/12.0,5.0/6.0},{1.0/12.0,1.0/12.0,-1.0/6.0}}, 186 U[3][2] = {{1,0},{1,10},{1,-1}}, 187 V[2][2] = {{1,0},{0,1}}, 188 S[2] = {1,0}, 189 F[2] = {1,0}, 190 Fembed[2] = {1,1-gamma}, 191 Ferror[2] = {0,1}, 192 Serror[2] = {1,0}; 193 ierr = TSGLEERegister(TSGLEE23,p,s,r,gamma,&A[0][0],&B[0][0],&U[0][0],&V[0][0],S,F,NULL,Fembed,Ferror,Serror,0,NULL);CHKERRQ(ierr); 194 } 195 { 196 /* y-y~ form */ 197 const PetscInt 198 p = 2, 199 s = 4, 200 r = 2; 201 const PetscReal 202 gamma = 0, 203 A[4][4] = {{0,0,0,0},{0.75,0,0,0},{0.25,29.0/60.0,0,0},{-21.0/44.0,145.0/44.0,-20.0/11.0,0}}, 204 B[2][4] = {{109.0/275.0,58.0/75.0,-37.0/110.0,1.0/6.0},{3.0/11.0,0,75.0/88.0,-1.0/8.0}}, 205 U[4][2] = {{0,1},{75.0/58.0,-17.0/58.0},{0,1},{0,1}}, 206 V[2][2] = {{1,0},{0,1}}, 207 S[2] = {1,1}, 208 F[2] = {1,0}, 209 Fembed[2] = {0,1}, 210 Ferror[2] = {-1.0/(1.0-gamma),1.0/(1.0-gamma)}, 211 Serror[2] = {1.0-gamma,1.0}; 212 ierr = TSGLEERegister(TSGLEE24,p,s,r,gamma,&A[0][0],&B[0][0],&U[0][0],&V[0][0],S,F,NULL,Fembed,Ferror,Serror,0,NULL);CHKERRQ(ierr); 213 } 214 { 215 /* y-y~ form */ 216 const PetscInt 217 p = 2, 218 s = 5, 219 r = 2; 220 const PetscReal 221 gamma = 0, 222 A[5][5] = {{0,0,0,0,0}, 223 {-0.94079244066783383269,0,0,0,0}, 224 {0.64228187778301907108,0.10915356933958500042,0,0,0}, 225 {-0.51764297742287450812,0.74414270351096040738,-0.71404164927824538121,0,0}, 226 {-0.44696561556825969206,-0.76768425657590196518,0.20111608138142987881,0.93828186737840469796,0}}, 227 B[2][5] = {{-0.029309178948150356153,-0.49671981884013874923,0.34275801517650053274,0.32941112623949194988,0.85385985637229662276}, 228 {0.78133219686062535272,0.074238691892675897635,0.57957363498384957966,-0.24638502829674959968,-0.18875949544040123033}}, 229 U[5][2] = {{0.16911424754448327735,0.83088575245551672265}, 230 {0.53638465733199574340,0.46361534266800425660}, 231 {0.39901579167169582526,0.60098420832830417474}, 232 {0.87689005530618575480,0.12310994469381424520}, 233 {0.99056100455550913009,0.0094389954444908699092}}, 234 V[2][2] = {{1,0},{0,1}}, 235 S[2] = {1,1}, 236 F[2] = {1,0}, 237 Fembed[2] = {0,1}, 238 Ferror[2] = {-1.0/(1.0-gamma),1.0/(1.0-gamma)}, 239 Serror[2] = {1.0-gamma,1.0}; 240 ierr = TSGLEERegister(TSGLEE25I,p,s,r,gamma,&A[0][0],&B[0][0],&U[0][0],&V[0][0],S,F,NULL,Fembed,Ferror,Serror,0,NULL);CHKERRQ(ierr); 241 } 242 { 243 /* y-y~ form */ 244 const PetscInt 245 p = 3, 246 s = 5, 247 r = 2; 248 const PetscReal 249 gamma = 0, 250 A[5][5] = {{0,0,0,0,0}, 251 {- 2169604947363702313.0 / 24313474998937147335.0,0,0,0,0}, 252 {46526746497697123895.0 / 94116917485856474137.0,-10297879244026594958.0 / 49199457603717988219.0,0,0,0}, 253 {23364788935845982499.0 / 87425311444725389446.0,-79205144337496116638.0 / 148994349441340815519.0,40051189859317443782.0 / 36487615018004984309.0,0,0}, 254 {42089522664062539205.0 / 124911313006412840286.0,-15074384760342762939.0 / 137927286865289746282.0,-62274678522253371016.0 / 125918573676298591413.0,13755475729852471739.0 / 79257927066651693390.0}}, 255 B[2][5] = {{61546696837458703723.0 / 56982519523786160813.0,-55810892792806293355.0 / 206957624151308356511.0,24061048952676379087.0 / 158739347956038723465.0,3577972206874351339.0 / 7599733370677197135.0,-59449832954780563947.0 / 137360038685338563670.0}, 256 {- 9738262186984159168.0 / 99299082461487742983.0,-32797097931948613195.0 / 61521565616362163366.0,42895514606418420631.0 / 71714201188501437336.0,22608567633166065068.0 / 55371917805607957003.0,94655809487476459565.0 / 151517167160302729021.0}}, 257 U[5][2] = {{70820309139834661559.0 / 80863923579509469826.0,10043614439674808267.0 / 80863923579509469826.0}, 258 {161694774978034105510.0 / 106187653640211060371.0,-55507121337823045139.0 / 106187653640211060371.0}, 259 {78486094644566264568.0 / 88171030896733822981.0,9684936252167558413.0 / 88171030896733822981.0}, 260 {65394922146334854435.0 / 84570853840405479554.0,19175931694070625119.0 / 84570853840405479554.0}, 261 {8607282770183754108.0 / 108658046436496925911.0,100050763666313171803.0 / 108658046436496925911.0}}, 262 V[2][2] = {{1,0},{0,1}}, 263 S[2] = {1,1}, 264 F[2] = {1,0}, 265 Fembed[2] = {0,1}, 266 Ferror[2] = {-1.0/(1.0-gamma),1.0/(1.0-gamma)}, 267 Serror[2] = {1.0-gamma,1.0}; 268 ierr = TSGLEERegister(TSGLEE35,p,s,r,gamma,&A[0][0],&B[0][0],&U[0][0],&V[0][0],S,F,NULL,Fembed,Ferror,Serror,0,NULL);CHKERRQ(ierr); 269 } 270 { 271 /* y-eps form */ 272 const PetscInt 273 p = 2, 274 s = 6, 275 r = 2; 276 const PetscReal 277 gamma = 0.25, 278 A[6][6] = {{0,0,0,0,0,0}, 279 {1,0,0,0,0,0}, 280 {0,0,0,0,0,0}, 281 {0,0,0.5,0,0,0}, 282 {0,0,0.25,0.25,0,0}, 283 {0,0,0.25,0.25,0.5,0}}, 284 B[2][6] = {{0.5,0.5,0,0,0,0}, 285 {-2.0/3.0,-2.0/3.0,1.0/3.0,1.0/3.0,1.0/3.0,1.0/3.0}}, 286 U[6][2] = {{1,0},{1,0},{1,0.75},{1,0.75},{1,0.75},{1,0.75}}, 287 V[2][2] = {{1,0},{0,1}}, 288 S[2] = {1,0}, 289 F[2] = {1,0}, 290 Fembed[2] = {1,1-gamma}, 291 Ferror[2] = {0,1}, 292 Serror[2] = {1,0}; 293 ierr = TSGLEERegister(TSGLEEEXRK2A,p,s,r,gamma,&A[0][0],&B[0][0],&U[0][0],&V[0][0],S,F,NULL,Fembed,Ferror,Serror,0,NULL);CHKERRQ(ierr); 294 } 295 { 296 /* y-eps form */ 297 const PetscInt 298 p = 3, 299 s = 8, 300 r = 2; 301 const PetscReal 302 gamma = 0, 303 A[8][8] = {{0,0,0,0,0,0,0,0}, 304 {0.5,0,0,0,0,0,0,0}, 305 {-1,2,0,0,0,0,0,0}, 306 {1.0/6.0,2.0/3.0,1.0/6.0,0,0,0,0,0}, 307 {0,0,0,0,0,0,0,0}, 308 {-7.0/24.0,1.0/3.0,1.0/12.0,-1.0/8.0,0.5,0,0,0}, 309 {7.0/6.0,-4.0/3.0,-1.0/3.0,0.5,-1.0,2.0,0,0}, 310 {0,0,0,0,1.0/6.0,2.0/3.0,1.0/6.0,0}}, 311 B[2][8] = {{1.0/6.0,2.0/3.0,1.0/6.0,0,0,0,0,0}, 312 {-1.0/6.0,-2.0/3.0,-1.0/6.0,0,1.0/6.0,2.0/3.0,1.0/6.0,0}}, 313 U[8][2] = {{1,0},{1,0},{1,0},{1,0},{1,1},{1,1},{1,1},{1,1}}, 314 V[2][2] = {{1,0},{0,1}}, 315 S[2] = {1,0}, 316 F[2] = {1,0}, 317 Fembed[2] = {1,1-gamma}, 318 Ferror[2] = {0,1}, 319 Serror[2] = {1,0}; 320 ierr = TSGLEERegister(TSGLEERK32G1,p,s,r,gamma,&A[0][0],&B[0][0],&U[0][0],&V[0][0],S,F,NULL,Fembed,Ferror,Serror,0,NULL);CHKERRQ(ierr); 321 } 322 { 323 /* y-eps form */ 324 const PetscInt 325 p = 2, 326 s = 9, 327 r = 2; 328 const PetscReal 329 gamma = 0.25, 330 A[9][9] = {{0,0,0,0,0,0,0,0,0}, 331 {0.585786437626904966,0,0,0,0,0,0,0,0}, 332 {0.149999999999999994,0.849999999999999978,0,0,0,0,0,0,0}, 333 {0,0,0,0,0,0,0,0,0}, 334 {0,0,0,0.292893218813452483,0,0,0,0,0}, 335 {0,0,0,0.0749999999999999972,0.424999999999999989,0,0,0,0}, 336 {0,0,0,0.176776695296636893,0.176776695296636893,0.146446609406726241,0,0,0}, 337 {0,0,0,0.176776695296636893,0.176776695296636893,0.146446609406726241,0.292893218813452483,0,0}, 338 {0,0,0,0.176776695296636893,0.176776695296636893,0.146446609406726241,0.0749999999999999972,0.424999999999999989,0}}, 339 B[2][9] = {{0.353553390593273786,0.353553390593273786,0.292893218813452483,0,0,0,0,0,0}, 340 {-0.471404520791031678,-0.471404520791031678,-0.390524291751269959,0.235702260395515839,0.235702260395515839,0.195262145875634979,0.235702260395515839,0.235702260395515839,0.195262145875634979}}, 341 U[9][2] = {{1,0},{1,0},{1,0},{1,0.75},{1,0.75},{1,0.75},{1,0.75},{1,0.75},{1,0.75}}, 342 V[2][2] = {{1,0},{0,1}}, 343 S[2] = {1,0}, 344 F[2] = {1,0}, 345 Fembed[2] = {1,1-gamma}, 346 Ferror[2] = {0,1}, 347 Serror[2] = {1,0}; 348 ierr = TSGLEERegister(TSGLEERK285EX,p,s,r,gamma,&A[0][0],&B[0][0],&U[0][0],&V[0][0],S,F,NULL,Fembed,Ferror,Serror,0,NULL);CHKERRQ(ierr); 349 } 350 351 PetscFunctionReturn(0); 352 } 353 354 #undef __FUNCT__ 355 #define __FUNCT__ "TSGLEERegisterDestroy" 356 /*@C 357 TSGLEERegisterDestroy - Frees the list of schemes that were registered by TSGLEERegister(). 358 359 Not Collective 360 361 Level: advanced 362 363 .keywords: TSGLEE, register, destroy 364 .seealso: TSGLEERegister(), TSGLEERegisterAll() 365 @*/ 366 PetscErrorCode TSGLEERegisterDestroy(void) 367 { 368 PetscErrorCode ierr; 369 GLEETableauLink link; 370 371 PetscFunctionBegin; 372 while ((link = GLEETableauList)) { 373 GLEETableau t = &link->tab; 374 GLEETableauList = link->next; 375 ierr = PetscFree5(t->A,t->B,t->U,t->V,t->c); CHKERRQ(ierr); 376 ierr = PetscFree2(t->S,t->F); CHKERRQ(ierr); 377 ierr = PetscFree (t->Fembed); CHKERRQ(ierr); 378 ierr = PetscFree (t->Ferror); CHKERRQ(ierr); 379 ierr = PetscFree (t->Serror); CHKERRQ(ierr); 380 ierr = PetscFree (t->binterp); CHKERRQ(ierr); 381 ierr = PetscFree (t->name); CHKERRQ(ierr); 382 ierr = PetscFree (link); CHKERRQ(ierr); 383 } 384 TSGLEERegisterAllCalled = PETSC_FALSE; 385 PetscFunctionReturn(0); 386 } 387 388 #undef __FUNCT__ 389 #define __FUNCT__ "TSGLEEInitializePackage" 390 /*@C 391 TSGLEEInitializePackage - This function initializes everything in the TSGLEE package. It is called 392 from PetscDLLibraryRegister() when using dynamic libraries, and on the first call to TSCreate_GLEE() 393 when using static libraries. 394 395 Level: developer 396 397 .keywords: TS, TSGLEE, initialize, package 398 .seealso: PetscInitialize() 399 @*/ 400 PetscErrorCode TSGLEEInitializePackage(void) 401 { 402 PetscErrorCode ierr; 403 404 PetscFunctionBegin; 405 if (TSGLEEPackageInitialized) PetscFunctionReturn(0); 406 TSGLEEPackageInitialized = PETSC_TRUE; 407 ierr = TSGLEERegisterAll();CHKERRQ(ierr); 408 ierr = PetscObjectComposedDataRegister(&explicit_stage_time_id);CHKERRQ(ierr); 409 ierr = PetscRegisterFinalize(TSGLEEFinalizePackage);CHKERRQ(ierr); 410 PetscFunctionReturn(0); 411 } 412 413 #undef __FUNCT__ 414 #define __FUNCT__ "TSGLEEFinalizePackage" 415 /*@C 416 TSGLEEFinalizePackage - This function destroys everything in the TSGLEE package. It is 417 called from PetscFinalize(). 418 419 Level: developer 420 421 .keywords: Petsc, destroy, package 422 .seealso: PetscFinalize() 423 @*/ 424 PetscErrorCode TSGLEEFinalizePackage(void) 425 { 426 PetscErrorCode ierr; 427 428 PetscFunctionBegin; 429 TSGLEEPackageInitialized = PETSC_FALSE; 430 ierr = TSGLEERegisterDestroy();CHKERRQ(ierr); 431 PetscFunctionReturn(0); 432 } 433 434 #undef __FUNCT__ 435 #define __FUNCT__ "TSGLEERegister" 436 /*@C 437 TSGLEERegister - register an GLEE scheme by providing the entries in the Butcher tableau 438 439 Not Collective, but the same schemes should be registered on all processes on which they will be used 440 441 Input Parameters: 442 + name - identifier for method 443 . order - order of method 444 . s - number of stages 445 . r - number of steps 446 . gamma - LTE ratio 447 . A - stage coefficients (dimension s*s, row-major) 448 . B - step completion coefficients (dimension r*s, row-major) 449 . U - method coefficients (dimension s*r, row-major) 450 . V - method coefficients (dimension r*r, row-major) 451 . S - starting coefficients 452 . F - finishing coefficients 453 . c - abscissa (dimension s; NULL to use row sums of A) 454 . Fembed - step completion coefficients for embedded method 455 . Ferror - error computation coefficients 456 . Serror - error initialization coefficients 457 . pinterp - order of interpolation (0 if unavailable) 458 - binterp - array of interpolation coefficients (NULL if unavailable) 459 460 Notes: 461 Several GLEE methods are provided, this function is only needed to create new methods. 462 463 Level: advanced 464 465 .keywords: TS, register 466 467 .seealso: TSGLEE 468 @*/ 469 PetscErrorCode TSGLEERegister(TSGLEEType name,PetscInt order,PetscInt s, PetscInt r, 470 PetscReal gamma, 471 const PetscReal A[],const PetscReal B[], 472 const PetscReal U[],const PetscReal V[], 473 const PetscReal S[],const PetscReal F[], 474 const PetscReal c[], 475 const PetscReal Fembed[],const PetscReal Ferror[], 476 const PetscReal Serror[], 477 PetscInt pinterp, const PetscReal binterp[]) 478 { 479 PetscErrorCode ierr; 480 GLEETableauLink link; 481 GLEETableau t; 482 PetscInt i,j; 483 484 PetscFunctionBegin; 485 ierr = PetscMalloc(sizeof(*link),&link);CHKERRQ(ierr); 486 ierr = PetscMemzero(link,sizeof(*link));CHKERRQ(ierr); 487 t = &link->tab; 488 ierr = PetscStrallocpy(name,&t->name);CHKERRQ(ierr); 489 t->order = order; 490 t->s = s; 491 t->r = r; 492 t->gamma = gamma; 493 ierr = PetscMalloc5(s*s,&t->A,r*r,&t->V,s,&t->c,r*s,&t->B,s*r,&t->U); CHKERRQ(ierr); 494 ierr = PetscMalloc2(r,&t->S,r,&t->F); CHKERRQ(ierr); 495 ierr = PetscMemcpy(t->A,A,s*s*sizeof(A[0]));CHKERRQ(ierr); 496 ierr = PetscMemcpy(t->B,B,r*s*sizeof(B[0]));CHKERRQ(ierr); 497 ierr = PetscMemcpy(t->U,U,s*r*sizeof(U[0]));CHKERRQ(ierr); 498 ierr = PetscMemcpy(t->V,V,r*r*sizeof(V[0]));CHKERRQ(ierr); 499 ierr = PetscMemcpy(t->S,S,r *sizeof(S[0]));CHKERRQ(ierr); 500 ierr = PetscMemcpy(t->F,F,r *sizeof(F[0]));CHKERRQ(ierr); 501 if (c) { 502 ierr = PetscMemcpy(t->c,c,s*sizeof(c[0])); CHKERRQ(ierr); 503 } else { 504 for (i=0; i<s; i++) for (j=0,t->c[i]=0; j<s; j++) t->c[i] += A[i*s+j]; 505 } 506 ierr = PetscMalloc1(r,&t->Fembed);CHKERRQ(ierr); 507 ierr = PetscMalloc1(r,&t->Ferror);CHKERRQ(ierr); 508 ierr = PetscMalloc1(r,&t->Serror);CHKERRQ(ierr); 509 ierr = PetscMemcpy(t->Fembed,Fembed,r*sizeof(Fembed[0]));CHKERRQ(ierr); 510 ierr = PetscMemcpy(t->Ferror,Ferror,r*sizeof(Ferror[0]));CHKERRQ(ierr); 511 ierr = PetscMemcpy(t->Serror,Serror,r*sizeof(Serror[0]));CHKERRQ(ierr); 512 t->pinterp = pinterp; 513 ierr = PetscMalloc1(s*pinterp,&t->binterp);CHKERRQ(ierr); 514 ierr = PetscMemcpy(t->binterp,binterp,s*pinterp*sizeof(binterp[0]));CHKERRQ(ierr); 515 516 link->next = GLEETableauList; 517 GLEETableauList = link; 518 PetscFunctionReturn(0); 519 } 520 521 #undef __FUNCT__ 522 #define __FUNCT__ "TSEvaluateStep_GLEE" 523 static PetscErrorCode TSEvaluateStep_GLEE(TS ts,PetscInt order,Vec X,PetscBool *done) 524 { 525 TS_GLEE *glee = (TS_GLEE*) ts->data; 526 GLEETableau tab = glee->tableau; 527 PetscReal h, *B = tab->B, *V = tab->V, 528 *F = tab->F, 529 *Fembed = tab->Fembed; 530 PetscInt s = tab->s, r = tab->r, i, j; 531 Vec *Y = glee->Y, *YdotStage = glee->YdotStage; 532 PetscScalar *w = glee->work; 533 PetscErrorCode ierr; 534 535 PetscFunctionBegin; 536 537 switch (glee->status) { 538 case TS_STEP_INCOMPLETE: 539 case TS_STEP_PENDING: 540 h = ts->time_step; break; 541 case TS_STEP_COMPLETE: 542 h = ts->ptime - ts->ptime_prev; break; 543 default: SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_PLIB,"Invalid TSStepStatus"); 544 } 545 546 if (order == tab->order) { 547 548 /* Note: Irrespective of whether status is TS_STEP_INCOMPLETE 549 or TS_STEP_COMPLETE, glee->X has the solution at the 550 beginning of the time step. So no need to roll-back. 551 */ 552 if (glee->status == TS_STEP_INCOMPLETE) { 553 for (i=0; i<r; i++) { 554 ierr = VecZeroEntries(Y[i]); CHKERRQ(ierr); 555 ierr = VecMAXPY(Y[i],r,(V+i*r),glee->X); CHKERRQ(ierr); 556 for (j=0; j<s; j++) w[j] = h*B[i*s+j]; 557 ierr = VecMAXPY(Y[i],s,w,YdotStage); CHKERRQ(ierr); 558 } 559 ierr = VecZeroEntries(X);CHKERRQ(ierr); 560 ierr = VecMAXPY(X,r,F,Y);CHKERRQ(ierr); 561 } else {ierr = VecCopy(ts->vec_sol,X);CHKERRQ(ierr);} 562 PetscFunctionReturn(0); 563 564 } else if (order == tab->order-1) { 565 566 /* Complete with the embedded method (Fembed) */ 567 for (i=0; i<r; i++) { 568 ierr = VecZeroEntries(Y[i]); CHKERRQ(ierr); 569 ierr = VecMAXPY(Y[i],r,(V+i*r),glee->X); CHKERRQ(ierr); 570 for (j=0; j<s; j++) w[j] = h*B[i*s+j]; 571 ierr = VecMAXPY(Y[i],s,w,YdotStage); CHKERRQ(ierr); 572 } 573 ierr = VecZeroEntries(X);CHKERRQ(ierr); 574 ierr = VecMAXPY(X,r,Fembed,Y);CHKERRQ(ierr); 575 576 if (done) *done = PETSC_TRUE; 577 PetscFunctionReturn(0); 578 } 579 if (done) *done = PETSC_FALSE; 580 else SETERRQ3(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"GLEE '%s' of order %D cannot evaluate step at order %D",tab->name,tab->order,order); 581 PetscFunctionReturn(0); 582 } 583 584 #undef __FUNCT__ 585 #define __FUNCT__ "TSStep_GLEE" 586 static PetscErrorCode TSStep_GLEE(TS ts) 587 { 588 TS_GLEE *glee = (TS_GLEE*)ts->data; 589 GLEETableau tab = glee->tableau; 590 const PetscInt s = tab->s, r = tab->r; 591 PetscReal *A = tab->A, *U = tab->U, 592 *F = tab->F, 593 *c = tab->c; 594 Vec *Y = glee->Y, *X = glee->X, 595 *YStage = glee->YStage, 596 *YdotStage = glee->YdotStage, 597 W = glee->W; 598 SNES snes; 599 PetscScalar *w = glee->work; 600 TSAdapt adapt; 601 PetscInt i,j,reject,next_scheme,its,lits; 602 PetscReal next_time_step; 603 PetscReal t; 604 PetscBool accept; 605 PetscErrorCode ierr; 606 607 PetscFunctionBegin; 608 ierr = PetscCitationsRegister(citation,&cited);CHKERRQ(ierr); 609 610 for (i=0; i<r; i++) { ierr = VecCopy(Y[i],X[i]); CHKERRQ(ierr); } 611 612 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 613 next_time_step = ts->time_step; 614 t = ts->ptime; 615 accept = PETSC_TRUE; 616 glee->status = TS_STEP_INCOMPLETE; 617 618 for (reject=0; reject<ts->max_reject && !ts->reason; reject++,ts->reject++) { 619 620 PetscReal h = ts->time_step; 621 ierr = TSPreStep(ts);CHKERRQ(ierr); 622 623 for (i=0; i<s; i++) { 624 625 glee->stage_time = t + h*c[i]; 626 ierr = TSPreStage(ts,glee->stage_time); CHKERRQ(ierr); 627 628 if (A[i*s+i] == 0) { /* Explicit stage */ 629 ierr = VecZeroEntries(YStage[i]);CHKERRQ(ierr); 630 ierr = VecMAXPY(YStage[i],r,(U+i*r),X);CHKERRQ(ierr); 631 for (j=0; j<i; j++) w[j] = h*A[i*s+j]; 632 ierr = VecMAXPY(YStage[i],i,w,YdotStage);CHKERRQ(ierr); 633 } else { /* Implicit stage */ 634 glee->scoeff = 1.0/A[i*s+i]; 635 /* compute right-hand-side */ 636 ierr = VecZeroEntries(W);CHKERRQ(ierr); 637 ierr = VecMAXPY(W,r,(U+i*r),X);CHKERRQ(ierr); 638 for (j=0; j<i; j++) w[j] = h*A[i*s+j]; 639 ierr = VecMAXPY(W,i,w,YdotStage);CHKERRQ(ierr); 640 ierr = VecScale(W,glee->scoeff/h);CHKERRQ(ierr); 641 /* set initial guess */ 642 ierr = VecCopy(i>0 ? YStage[i-1] : ts->vec_sol,YStage[i]);CHKERRQ(ierr); 643 /* solve for this stage */ 644 ierr = SNESSolve(snes,W,YStage[i]);CHKERRQ(ierr); 645 ierr = SNESGetIterationNumber(snes,&its);CHKERRQ(ierr); 646 ierr = SNESGetLinearSolveIterations(snes,&lits);CHKERRQ(ierr); 647 ts->snes_its += its; ts->ksp_its += lits; 648 } 649 ierr = TSGetAdapt(ts,&adapt);CHKERRQ(ierr); 650 ierr = TSAdaptCheckStage(adapt,ts,glee->stage_time,YStage[i],&accept);CHKERRQ(ierr); 651 if (!accept) goto reject_step; 652 ierr = TSPostStage(ts,glee->stage_time,i,YStage); CHKERRQ(ierr); 653 ierr = TSComputeRHSFunction(ts,t+h*c[i],YStage[i],YdotStage[i]);CHKERRQ(ierr); 654 } 655 ierr = TSEvaluateStep(ts,tab->order,ts->vec_sol,NULL);CHKERRQ(ierr); 656 glee->status = TS_STEP_PENDING; 657 658 /* Register only the current method as a candidate because we're not supporting multiple candidates yet. */ 659 ierr = TSGetAdapt(ts,&adapt);CHKERRQ(ierr); 660 ierr = TSAdaptCandidatesClear(adapt);CHKERRQ(ierr); 661 ierr = TSAdaptCandidateAdd(adapt,tab->name,tab->order,1,tab->ccfl,1.*tab->s,PETSC_TRUE);CHKERRQ(ierr); 662 ierr = TSAdaptChoose(adapt,ts,ts->time_step,&next_scheme,&next_time_step,&accept);CHKERRQ(ierr); 663 if (accept) { 664 /* ignore next_scheme for now */ 665 ts->ptime += ts->time_step; 666 ts->time_step = next_time_step; 667 glee->status = TS_STEP_COMPLETE; 668 /* compute and store the global error */ 669 /* Note: this is not needed if TSAdaptGLEE is not used */ 670 ierr = TSGetTimeError(ts,0,&(glee->yGErr));CHKERRQ(ierr); 671 ierr = PetscObjectComposedDataSetReal((PetscObject)ts->vec_sol,explicit_stage_time_id,ts->ptime);CHKERRQ(ierr); 672 break; 673 } else { /* Roll back the current step */ 674 ierr = VecMAXPY(ts->vec_sol,r,F,X);CHKERRQ(ierr); 675 ts->time_step = next_time_step; 676 glee->status = TS_STEP_INCOMPLETE; 677 } 678 reject_step: continue; 679 } 680 if (glee->status != TS_STEP_COMPLETE && !ts->reason) ts->reason = TS_DIVERGED_STEP_REJECTED; 681 PetscFunctionReturn(0); 682 } 683 684 #undef __FUNCT__ 685 #define __FUNCT__ "TSInterpolate_GLEE" 686 static PetscErrorCode TSInterpolate_GLEE(TS ts,PetscReal itime,Vec X) 687 { 688 TS_GLEE *glee = (TS_GLEE*)ts->data; 689 PetscInt s=glee->tableau->s, pinterp=glee->tableau->pinterp,i,j; 690 PetscReal h,tt,t; 691 PetscScalar *b; 692 const PetscReal *B = glee->tableau->binterp; 693 PetscErrorCode ierr; 694 695 PetscFunctionBegin; 696 if (!B) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"TSGLEE %s does not have an interpolation formula",glee->tableau->name); 697 switch (glee->status) { 698 case TS_STEP_INCOMPLETE: 699 case TS_STEP_PENDING: 700 h = ts->time_step; 701 t = (itime - ts->ptime)/h; 702 break; 703 case TS_STEP_COMPLETE: 704 h = ts->ptime - ts->ptime_prev; 705 t = (itime - ts->ptime)/h + 1; /* In the interval [0,1] */ 706 break; 707 default: SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_PLIB,"Invalid TSStepStatus"); 708 } 709 ierr = PetscMalloc1(s,&b);CHKERRQ(ierr); 710 for (i=0; i<s; i++) b[i] = 0; 711 for (j=0,tt=t; j<pinterp; j++,tt*=t) { 712 for (i=0; i<s; i++) { 713 b[i] += h * B[i*pinterp+j] * tt; 714 } 715 } 716 ierr = VecCopy(glee->YStage[0],X);CHKERRQ(ierr); 717 ierr = VecMAXPY(X,s,b,glee->YdotStage);CHKERRQ(ierr); 718 ierr = PetscFree(b);CHKERRQ(ierr); 719 PetscFunctionReturn(0); 720 } 721 722 /*------------------------------------------------------------*/ 723 #undef __FUNCT__ 724 #define __FUNCT__ "TSReset_GLEE" 725 static PetscErrorCode TSReset_GLEE(TS ts) 726 { 727 TS_GLEE *glee = (TS_GLEE*)ts->data; 728 PetscInt s, r; 729 PetscErrorCode ierr; 730 731 PetscFunctionBegin; 732 if (!glee->tableau) PetscFunctionReturn(0); 733 s = glee->tableau->s; 734 r = glee->tableau->r; 735 ierr = VecDestroyVecs(r,&glee->Y);CHKERRQ(ierr); 736 ierr = VecDestroyVecs(r,&glee->X);CHKERRQ(ierr); 737 ierr = VecDestroyVecs(s,&glee->YStage);CHKERRQ(ierr); 738 ierr = VecDestroyVecs(s,&glee->YdotStage);CHKERRQ(ierr); 739 ierr = VecDestroy(&glee->Ydot);CHKERRQ(ierr); 740 ierr = VecDestroy(&glee->yGErr);CHKERRQ(ierr); 741 ierr = VecDestroy(&glee->W);CHKERRQ(ierr); 742 ierr = PetscFree(glee->work);CHKERRQ(ierr); 743 PetscFunctionReturn(0); 744 } 745 746 #undef __FUNCT__ 747 #define __FUNCT__ "TSDestroy_GLEE" 748 static PetscErrorCode TSDestroy_GLEE(TS ts) 749 { 750 PetscErrorCode ierr; 751 752 PetscFunctionBegin; 753 ierr = TSReset_GLEE(ts);CHKERRQ(ierr); 754 ierr = PetscFree(ts->data);CHKERRQ(ierr); 755 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSGLEEGetType_C",NULL);CHKERRQ(ierr); 756 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSGLEESetType_C",NULL);CHKERRQ(ierr); 757 PetscFunctionReturn(0); 758 } 759 760 #undef __FUNCT__ 761 #define __FUNCT__ "TSGLEEGetVecs" 762 static PetscErrorCode TSGLEEGetVecs(TS ts,DM dm,Vec *Ydot) 763 { 764 TS_GLEE *glee = (TS_GLEE*)ts->data; 765 PetscErrorCode ierr; 766 767 PetscFunctionBegin; 768 if (Ydot) { 769 if (dm && dm != ts->dm) { 770 ierr = DMGetNamedGlobalVector(dm,"TSGLEE_Ydot",Ydot);CHKERRQ(ierr); 771 } else *Ydot = glee->Ydot; 772 } 773 PetscFunctionReturn(0); 774 } 775 776 777 #undef __FUNCT__ 778 #define __FUNCT__ "TSGLEERestoreVecs" 779 static PetscErrorCode TSGLEERestoreVecs(TS ts,DM dm,Vec *Ydot) 780 { 781 PetscErrorCode ierr; 782 783 PetscFunctionBegin; 784 if (Ydot) { 785 if (dm && dm != ts->dm) { 786 ierr = DMRestoreNamedGlobalVector(dm,"TSGLEE_Ydot",Ydot);CHKERRQ(ierr); 787 } 788 } 789 PetscFunctionReturn(0); 790 } 791 792 /* 793 This defines the nonlinear equation that is to be solved with SNES 794 */ 795 #undef __FUNCT__ 796 #define __FUNCT__ "SNESTSFormFunction_GLEE" 797 static PetscErrorCode SNESTSFormFunction_GLEE(SNES snes,Vec X,Vec F,TS ts) 798 { 799 TS_GLEE *glee = (TS_GLEE*)ts->data; 800 DM dm,dmsave; 801 Vec Ydot; 802 PetscReal shift = glee->scoeff / ts->time_step; 803 PetscErrorCode ierr; 804 805 PetscFunctionBegin; 806 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 807 ierr = TSGLEEGetVecs(ts,dm,&Ydot);CHKERRQ(ierr); 808 /* Set Ydot = shift*X */ 809 ierr = VecCopy(X,Ydot);CHKERRQ(ierr); 810 ierr = VecScale(Ydot,shift);CHKERRQ(ierr); 811 dmsave = ts->dm; 812 ts->dm = dm; 813 814 ierr = TSComputeIFunction(ts,glee->stage_time,X,Ydot,F,PETSC_FALSE);CHKERRQ(ierr); 815 816 ts->dm = dmsave; 817 ierr = TSGLEERestoreVecs(ts,dm,&Ydot);CHKERRQ(ierr); 818 PetscFunctionReturn(0); 819 } 820 821 #undef __FUNCT__ 822 #define __FUNCT__ "SNESTSFormJacobian_GLEE" 823 static PetscErrorCode SNESTSFormJacobian_GLEE(SNES snes,Vec X,Mat A,Mat B,TS ts) 824 { 825 TS_GLEE *glee = (TS_GLEE*)ts->data; 826 DM dm,dmsave; 827 Vec Ydot; 828 PetscReal shift = glee->scoeff / ts->time_step; 829 PetscErrorCode ierr; 830 831 PetscFunctionBegin; 832 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 833 ierr = TSGLEEGetVecs(ts,dm,&Ydot);CHKERRQ(ierr); 834 /* glee->Ydot has already been computed in SNESTSFormFunction_GLEE (SNES guarantees this) */ 835 dmsave = ts->dm; 836 ts->dm = dm; 837 838 ierr = TSComputeIJacobian(ts,glee->stage_time,X,Ydot,shift,A,B,PETSC_FALSE);CHKERRQ(ierr); 839 840 ts->dm = dmsave; 841 ierr = TSGLEERestoreVecs(ts,dm,&Ydot);CHKERRQ(ierr); 842 PetscFunctionReturn(0); 843 } 844 845 #undef __FUNCT__ 846 #define __FUNCT__ "DMCoarsenHook_TSGLEE" 847 static PetscErrorCode DMCoarsenHook_TSGLEE(DM fine,DM coarse,void *ctx) 848 { 849 PetscFunctionBegin; 850 PetscFunctionReturn(0); 851 } 852 853 #undef __FUNCT__ 854 #define __FUNCT__ "DMRestrictHook_TSGLEE" 855 static PetscErrorCode DMRestrictHook_TSGLEE(DM fine,Mat restrct,Vec rscale,Mat inject,DM coarse,void *ctx) 856 { 857 PetscFunctionBegin; 858 PetscFunctionReturn(0); 859 } 860 861 862 #undef __FUNCT__ 863 #define __FUNCT__ "DMSubDomainHook_TSGLEE" 864 static PetscErrorCode DMSubDomainHook_TSGLEE(DM dm,DM subdm,void *ctx) 865 { 866 PetscFunctionBegin; 867 PetscFunctionReturn(0); 868 } 869 870 #undef __FUNCT__ 871 #define __FUNCT__ "DMSubDomainRestrictHook_TSGLEE" 872 static PetscErrorCode DMSubDomainRestrictHook_TSGLEE(DM dm,VecScatter gscat,VecScatter lscat,DM subdm,void *ctx) 873 { 874 PetscFunctionBegin; 875 PetscFunctionReturn(0); 876 } 877 878 #undef __FUNCT__ 879 #define __FUNCT__ "TSSetUp_GLEE" 880 static PetscErrorCode TSSetUp_GLEE(TS ts) 881 { 882 TS_GLEE *glee = (TS_GLEE*)ts->data; 883 GLEETableau tab; 884 PetscInt s,r; 885 PetscErrorCode ierr; 886 DM dm; 887 888 PetscFunctionBegin; 889 if (!glee->tableau) { 890 ierr = TSGLEESetType(ts,TSGLEEDefaultType);CHKERRQ(ierr); 891 } 892 tab = glee->tableau; 893 s = tab->s; 894 r = tab->r; 895 ierr = VecDuplicateVecs(ts->vec_sol,r,&glee->Y);CHKERRQ(ierr); 896 ierr = VecDuplicateVecs(ts->vec_sol,r,&glee->X);CHKERRQ(ierr); 897 ierr = VecDuplicateVecs(ts->vec_sol,s,&glee->YStage);CHKERRQ(ierr); 898 ierr = VecDuplicateVecs(ts->vec_sol,s,&glee->YdotStage);CHKERRQ(ierr); 899 ierr = VecDuplicate(ts->vec_sol,&glee->Ydot);CHKERRQ(ierr); 900 ierr = VecDuplicate(ts->vec_sol,&glee->yGErr);CHKERRQ(ierr); 901 ierr = VecZeroEntries(glee->yGErr);CHKERRQ(ierr); 902 ierr = VecDuplicate(ts->vec_sol,&glee->W);CHKERRQ(ierr); 903 ierr = PetscMalloc1(s,&glee->work);CHKERRQ(ierr); 904 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 905 if (dm) { 906 ierr = DMCoarsenHookAdd(dm,DMCoarsenHook_TSGLEE,DMRestrictHook_TSGLEE,ts);CHKERRQ(ierr); 907 ierr = DMSubDomainHookAdd(dm,DMSubDomainHook_TSGLEE,DMSubDomainRestrictHook_TSGLEE,ts);CHKERRQ(ierr); 908 } 909 PetscFunctionReturn(0); 910 } 911 912 #undef __FUNCT__ 913 #define __FUNCT__ "TSStartingMethod_GLEE" 914 PetscErrorCode TSStartingMethod_GLEE(TS ts) 915 { 916 TS_GLEE *glee = (TS_GLEE*)ts->data; 917 GLEETableau tab = glee->tableau; 918 PetscInt r=tab->r,i; 919 PetscReal *S=tab->S; 920 PetscErrorCode ierr; 921 922 PetscFunctionBegin; 923 for (i=0; i<r; i++) { 924 ierr = VecZeroEntries(glee->Y[i]);CHKERRQ(ierr); 925 ierr = VecAXPY(glee->Y[i],S[i],ts->vec_sol);CHKERRQ(ierr); 926 } 927 928 PetscFunctionReturn(0); 929 } 930 931 932 /*------------------------------------------------------------*/ 933 934 #undef __FUNCT__ 935 #define __FUNCT__ "TSSetFromOptions_GLEE" 936 static PetscErrorCode TSSetFromOptions_GLEE(PetscOptionItems *PetscOptionsObject,TS ts) 937 { 938 PetscErrorCode ierr; 939 char gleetype[256]; 940 941 PetscFunctionBegin; 942 ierr = PetscOptionsHead(PetscOptionsObject,"GLEE ODE solver options");CHKERRQ(ierr); 943 { 944 GLEETableauLink link; 945 PetscInt count,choice; 946 PetscBool flg; 947 const char **namelist; 948 949 ierr = PetscStrncpy(gleetype,TSGLEEDefaultType,sizeof(gleetype));CHKERRQ(ierr); 950 for (link=GLEETableauList,count=0; link; link=link->next,count++) ; 951 ierr = PetscMalloc1(count,&namelist);CHKERRQ(ierr); 952 for (link=GLEETableauList,count=0; link; link=link->next,count++) namelist[count] = link->tab.name; 953 ierr = PetscOptionsEList("-ts_glee_type","Family of GLEE method","TSGLEESetType",(const char*const*)namelist,count,gleetype,&choice,&flg);CHKERRQ(ierr); 954 ierr = TSGLEESetType(ts,flg ? namelist[choice] : gleetype);CHKERRQ(ierr); 955 ierr = PetscFree(namelist);CHKERRQ(ierr); 956 } 957 ierr = PetscOptionsTail();CHKERRQ(ierr); 958 PetscFunctionReturn(0); 959 } 960 961 #undef __FUNCT__ 962 #define __FUNCT__ "PetscFormatRealArray" 963 static PetscErrorCode PetscFormatRealArray(char buf[],size_t len,const char *fmt,PetscInt n,const PetscReal x[]) 964 { 965 PetscErrorCode ierr; 966 PetscInt i; 967 size_t left,count; 968 char *p; 969 970 PetscFunctionBegin; 971 for (i=0,p=buf,left=len; i<n; i++) { 972 ierr = PetscSNPrintfCount(p,left,fmt,&count,x[i]);CHKERRQ(ierr); 973 if (count >= left) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Insufficient space in buffer"); 974 left -= count; 975 p += count; 976 *p++ = ' '; 977 } 978 p[i ? 0 : -1] = 0; 979 PetscFunctionReturn(0); 980 } 981 982 #undef __FUNCT__ 983 #define __FUNCT__ "TSView_GLEE" 984 static PetscErrorCode TSView_GLEE(TS ts,PetscViewer viewer) 985 { 986 TS_GLEE *glee = (TS_GLEE*)ts->data; 987 GLEETableau tab = glee->tableau; 988 PetscBool iascii; 989 PetscErrorCode ierr; 990 TSAdapt adapt; 991 992 PetscFunctionBegin; 993 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 994 if (iascii) { 995 TSGLEEType gleetype; 996 char buf[512]; 997 ierr = TSGLEEGetType(ts,&gleetype);CHKERRQ(ierr); 998 ierr = PetscViewerASCIIPrintf(viewer," GLEE %s\n",gleetype);CHKERRQ(ierr); 999 ierr = PetscFormatRealArray(buf,sizeof(buf),"% 8.6f",tab->s,tab->c);CHKERRQ(ierr); 1000 ierr = PetscViewerASCIIPrintf(viewer," Abscissa c = %s\n",buf);CHKERRQ(ierr); 1001 /* Note: print out r as well */ 1002 } 1003 ierr = TSGetAdapt(ts,&adapt);CHKERRQ(ierr); 1004 ierr = TSAdaptView(adapt,viewer);CHKERRQ(ierr); 1005 PetscFunctionReturn(0); 1006 } 1007 1008 #undef __FUNCT__ 1009 #define __FUNCT__ "TSLoad_GLEE" 1010 static PetscErrorCode TSLoad_GLEE(TS ts,PetscViewer viewer) 1011 { 1012 PetscErrorCode ierr; 1013 SNES snes; 1014 TSAdapt tsadapt; 1015 1016 PetscFunctionBegin; 1017 ierr = TSGetAdapt(ts,&tsadapt);CHKERRQ(ierr); 1018 ierr = TSAdaptLoad(tsadapt,viewer);CHKERRQ(ierr); 1019 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 1020 ierr = SNESLoad(snes,viewer);CHKERRQ(ierr); 1021 /* function and Jacobian context for SNES when used with TS is always ts object */ 1022 ierr = SNESSetFunction(snes,NULL,NULL,ts);CHKERRQ(ierr); 1023 ierr = SNESSetJacobian(snes,NULL,NULL,NULL,ts);CHKERRQ(ierr); 1024 PetscFunctionReturn(0); 1025 } 1026 1027 #undef __FUNCT__ 1028 #define __FUNCT__ "TSGLEESetType" 1029 /*@C 1030 TSGLEESetType - Set the type of GLEE scheme 1031 1032 Logically collective 1033 1034 Input Parameter: 1035 + ts - timestepping context 1036 - gleetype - type of GLEE-scheme 1037 1038 Level: intermediate 1039 1040 .seealso: TSGLEEGetType(), TSGLEE 1041 @*/ 1042 PetscErrorCode TSGLEESetType(TS ts,TSGLEEType gleetype) 1043 { 1044 PetscErrorCode ierr; 1045 1046 PetscFunctionBegin; 1047 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1048 ierr = PetscTryMethod(ts,"TSGLEESetType_C",(TS,TSGLEEType),(ts,gleetype));CHKERRQ(ierr); 1049 PetscFunctionReturn(0); 1050 } 1051 1052 #undef __FUNCT__ 1053 #define __FUNCT__ "TSGLEEGetType" 1054 /*@C 1055 TSGLEEGetType - Get the type of GLEE scheme 1056 1057 Logically collective 1058 1059 Input Parameter: 1060 . ts - timestepping context 1061 1062 Output Parameter: 1063 . gleetype - type of GLEE-scheme 1064 1065 Level: intermediate 1066 1067 .seealso: TSGLEESetType() 1068 @*/ 1069 PetscErrorCode TSGLEEGetType(TS ts,TSGLEEType *gleetype) 1070 { 1071 PetscErrorCode ierr; 1072 1073 PetscFunctionBegin; 1074 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1075 ierr = PetscUseMethod(ts,"TSGLEEGetType_C",(TS,TSGLEEType*),(ts,gleetype));CHKERRQ(ierr); 1076 PetscFunctionReturn(0); 1077 } 1078 1079 #undef __FUNCT__ 1080 #define __FUNCT__ "TSGLEEGetType_GLEE" 1081 PetscErrorCode TSGLEEGetType_GLEE(TS ts,TSGLEEType *gleetype) 1082 { 1083 TS_GLEE *glee = (TS_GLEE*)ts->data; 1084 PetscErrorCode ierr; 1085 1086 PetscFunctionBegin; 1087 if (!glee->tableau) { 1088 ierr = TSGLEESetType(ts,TSGLEEDefaultType);CHKERRQ(ierr); 1089 } 1090 *gleetype = glee->tableau->name; 1091 PetscFunctionReturn(0); 1092 } 1093 #undef __FUNCT__ 1094 #define __FUNCT__ "TSGLEESetType_GLEE" 1095 PetscErrorCode TSGLEESetType_GLEE(TS ts,TSGLEEType gleetype) 1096 { 1097 TS_GLEE *glee = (TS_GLEE*)ts->data; 1098 PetscErrorCode ierr; 1099 PetscBool match; 1100 GLEETableauLink link; 1101 1102 PetscFunctionBegin; 1103 if (glee->tableau) { 1104 ierr = PetscStrcmp(glee->tableau->name,gleetype,&match);CHKERRQ(ierr); 1105 if (match) PetscFunctionReturn(0); 1106 } 1107 for (link = GLEETableauList; link; link=link->next) { 1108 ierr = PetscStrcmp(link->tab.name,gleetype,&match);CHKERRQ(ierr); 1109 if (match) { 1110 ierr = TSReset_GLEE(ts);CHKERRQ(ierr); 1111 glee->tableau = &link->tab; 1112 PetscFunctionReturn(0); 1113 } 1114 } 1115 SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_UNKNOWN_TYPE,"Could not find '%s'",gleetype); 1116 PetscFunctionReturn(0); 1117 } 1118 1119 #undef __FUNCT__ 1120 #define __FUNCT__ "TSGetStages_GLEE" 1121 static PetscErrorCode TSGetStages_GLEE(TS ts,PetscInt *ns,Vec **Y) 1122 { 1123 TS_GLEE *glee = (TS_GLEE*)ts->data; 1124 1125 PetscFunctionBegin; 1126 *ns = glee->tableau->s; 1127 if(Y) *Y = glee->YStage; 1128 PetscFunctionReturn(0); 1129 } 1130 1131 #undef __FUNCT__ 1132 #define __FUNCT__ "TSGetSolutionComponents_GLEE" 1133 PetscErrorCode TSGetSolutionComponents_GLEE(TS ts,PetscInt *n,Vec *Y) 1134 { 1135 TS_GLEE *glee = (TS_GLEE*)ts->data; 1136 GLEETableau tab = glee->tableau; 1137 PetscErrorCode ierr; 1138 1139 PetscFunctionBegin; 1140 if (!Y) *n = tab->r; 1141 else { 1142 if ((*n >= 0) && (*n < tab->r)) { 1143 ierr = VecCopy(glee->Y[*n],*Y); CHKERRQ(ierr); 1144 } else SETERRQ3(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_OUTOFRANGE,"Second argument (%d) out of range[%d,%d].",*n,0,tab->r-1); 1145 } 1146 PetscFunctionReturn(0); 1147 } 1148 1149 #undef __FUNCT__ 1150 #define __FUNCT__ "TSGetAuxSolution_GLEE" 1151 PetscErrorCode TSGetAuxSolution_GLEE(TS ts,Vec *X) 1152 { 1153 TS_GLEE *glee = (TS_GLEE*)ts->data; 1154 GLEETableau tab = glee->tableau; 1155 PetscReal *F = tab->Fembed; 1156 PetscInt r = tab->r; 1157 Vec *Y = glee->Y; 1158 PetscErrorCode ierr; 1159 1160 PetscFunctionBegin; 1161 ierr = VecZeroEntries(*X);CHKERRQ(ierr); 1162 ierr = VecMAXPY((*X),r,F,Y);CHKERRQ(ierr); 1163 PetscFunctionReturn(0); 1164 } 1165 1166 #undef __FUNCT__ 1167 #define __FUNCT__ "TSGetTimeError_GLEE" 1168 PetscErrorCode TSGetTimeError_GLEE(TS ts,PetscInt n,Vec *X) 1169 { 1170 TS_GLEE *glee = (TS_GLEE*)ts->data; 1171 GLEETableau tab = glee->tableau; 1172 PetscReal *F = tab->Ferror; 1173 PetscInt r = tab->r; 1174 Vec *Y = glee->Y; 1175 PetscErrorCode ierr; 1176 1177 PetscFunctionBegin; 1178 ierr = VecZeroEntries(*X);CHKERRQ(ierr); 1179 if(n==0){ 1180 ierr = VecMAXPY((*X),r,F,Y);CHKERRQ(ierr); 1181 } else if(n==-1) { 1182 *X=glee->yGErr; 1183 } 1184 PetscFunctionReturn(0); 1185 } 1186 1187 #undef __FUNCT__ 1188 #define __FUNCT__ "TSSetTimeError_GLEE" 1189 PetscErrorCode TSSetTimeError_GLEE(TS ts,Vec X) 1190 { 1191 TS_GLEE *glee = (TS_GLEE*)ts->data; 1192 GLEETableau tab = glee->tableau; 1193 PetscReal *S = tab->Serror; 1194 PetscInt r = tab->r,i; 1195 Vec *Y = glee->Y; 1196 PetscErrorCode ierr; 1197 1198 PetscFunctionBegin; 1199 if (r != 2) SETERRQ2(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"TSSetTimeError_GLEE not supported for '%s' with r=%D.",tab->name,tab->r); 1200 for (i=1; i<r; i++) { 1201 ierr = VecCopy(ts->vec_sol,Y[i]); CHKERRQ(ierr); 1202 ierr = VecAXPBY(Y[i],S[0],S[1],X); CHKERRQ(ierr); 1203 ierr = VecCopy(X,glee->yGErr); CHKERRQ(ierr); 1204 } 1205 PetscFunctionReturn(0); 1206 } 1207 1208 /* ------------------------------------------------------------ */ 1209 /*MC 1210 TSGLEE - ODE and DAE solver using General Linear with Error Estimation schemes 1211 1212 The user should provide the right hand side of the equation 1213 using TSSetRHSFunction(). 1214 1215 Notes: 1216 The default is TSGLEE35, it can be changed with TSGLEESetType() or -ts_glee_type 1217 1218 Level: beginner 1219 1220 .seealso: TSCreate(), TS, TSSetType(), TSGLEESetType(), TSGLEEGetType(), 1221 TSGLEE23, TTSGLEE24, TSGLEE35, TSGLEE25I, TSGLEEEXRK2A, 1222 TSGLEERK32G1, TSGLEERK285EX, TSGLEEType, TSGLEERegister() 1223 1224 M*/ 1225 #undef __FUNCT__ 1226 #define __FUNCT__ "TSCreate_GLEE" 1227 PETSC_EXTERN PetscErrorCode TSCreate_GLEE(TS ts) 1228 { 1229 TS_GLEE *th; 1230 PetscErrorCode ierr; 1231 1232 PetscFunctionBegin; 1233 #if !defined(PETSC_USE_DYNAMIC_LIBRARIES) 1234 ierr = TSGLEEInitializePackage();CHKERRQ(ierr); 1235 #endif 1236 1237 ts->ops->reset = TSReset_GLEE; 1238 ts->ops->destroy = TSDestroy_GLEE; 1239 ts->ops->view = TSView_GLEE; 1240 ts->ops->load = TSLoad_GLEE; 1241 ts->ops->setup = TSSetUp_GLEE; 1242 ts->ops->step = TSStep_GLEE; 1243 ts->ops->interpolate = TSInterpolate_GLEE; 1244 ts->ops->evaluatestep = TSEvaluateStep_GLEE; 1245 ts->ops->setfromoptions = TSSetFromOptions_GLEE; 1246 ts->ops->getstages = TSGetStages_GLEE; 1247 ts->ops->snesfunction = SNESTSFormFunction_GLEE; 1248 ts->ops->snesjacobian = SNESTSFormJacobian_GLEE; 1249 ts->ops->getsolutioncomponents = TSGetSolutionComponents_GLEE; 1250 ts->ops->getauxsolution = TSGetAuxSolution_GLEE; 1251 ts->ops->gettimeerror = TSGetTimeError_GLEE; 1252 ts->ops->settimeerror = TSSetTimeError_GLEE; 1253 ts->ops->startingmethod = TSStartingMethod_GLEE; 1254 1255 ierr = PetscNewLog(ts,&th);CHKERRQ(ierr); 1256 ts->data = (void*)th; 1257 1258 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSGLEEGetType_C",TSGLEEGetType_GLEE);CHKERRQ(ierr); 1259 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSGLEESetType_C",TSGLEESetType_GLEE);CHKERRQ(ierr); 1260 PetscFunctionReturn(0); 1261 } 1262