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,Y[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 ts->steps++; 668 glee->status = TS_STEP_COMPLETE; 669 /* compute and store the global error */ 670 /* Note: this is not needed if TSAdaptGLEE is not used */ 671 ierr = TSGetTimeError(ts,0,&(glee->yGErr));CHKERRQ(ierr); 672 ierr = PetscObjectComposedDataSetReal((PetscObject)ts->vec_sol,explicit_stage_time_id,ts->ptime);CHKERRQ(ierr); 673 break; 674 } else { /* Roll back the current step */ 675 ierr = VecMAXPY(ts->vec_sol,r,F,X);CHKERRQ(ierr); 676 ts->time_step = next_time_step; 677 glee->status = TS_STEP_INCOMPLETE; 678 } 679 reject_step: continue; 680 } 681 if (glee->status != TS_STEP_COMPLETE && !ts->reason) ts->reason = TS_DIVERGED_STEP_REJECTED; 682 PetscFunctionReturn(0); 683 } 684 685 #undef __FUNCT__ 686 #define __FUNCT__ "TSInterpolate_GLEE" 687 static PetscErrorCode TSInterpolate_GLEE(TS ts,PetscReal itime,Vec X) 688 { 689 TS_GLEE *glee = (TS_GLEE*)ts->data; 690 PetscInt s=glee->tableau->s, pinterp=glee->tableau->pinterp,i,j; 691 PetscReal h,tt,t; 692 PetscScalar *b; 693 const PetscReal *B = glee->tableau->binterp; 694 PetscErrorCode ierr; 695 696 PetscFunctionBegin; 697 if (!B) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"TSGLEE %s does not have an interpolation formula",glee->tableau->name); 698 switch (glee->status) { 699 case TS_STEP_INCOMPLETE: 700 case TS_STEP_PENDING: 701 h = ts->time_step; 702 t = (itime - ts->ptime)/h; 703 break; 704 case TS_STEP_COMPLETE: 705 h = ts->ptime - ts->ptime_prev; 706 t = (itime - ts->ptime)/h + 1; /* In the interval [0,1] */ 707 break; 708 default: SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_PLIB,"Invalid TSStepStatus"); 709 } 710 ierr = PetscMalloc1(s,&b);CHKERRQ(ierr); 711 for (i=0; i<s; i++) b[i] = 0; 712 for (j=0,tt=t; j<pinterp; j++,tt*=t) { 713 for (i=0; i<s; i++) { 714 b[i] += h * B[i*pinterp+j] * tt; 715 } 716 } 717 ierr = VecCopy(glee->YStage[0],X);CHKERRQ(ierr); 718 ierr = VecMAXPY(X,s,b,glee->YdotStage);CHKERRQ(ierr); 719 ierr = PetscFree(b);CHKERRQ(ierr); 720 PetscFunctionReturn(0); 721 } 722 723 /*------------------------------------------------------------*/ 724 #undef __FUNCT__ 725 #define __FUNCT__ "TSReset_GLEE" 726 static PetscErrorCode TSReset_GLEE(TS ts) 727 { 728 TS_GLEE *glee = (TS_GLEE*)ts->data; 729 PetscInt s, r; 730 PetscErrorCode ierr; 731 732 PetscFunctionBegin; 733 if (!glee->tableau) PetscFunctionReturn(0); 734 s = glee->tableau->s; 735 r = glee->tableau->r; 736 ierr = VecDestroyVecs(r,&glee->Y);CHKERRQ(ierr); 737 ierr = VecDestroyVecs(r,&glee->X);CHKERRQ(ierr); 738 ierr = VecDestroyVecs(s,&glee->YStage);CHKERRQ(ierr); 739 ierr = VecDestroyVecs(s,&glee->YdotStage);CHKERRQ(ierr); 740 ierr = VecDestroy(&glee->Ydot);CHKERRQ(ierr); 741 ierr = VecDestroy(&glee->yGErr);CHKERRQ(ierr); 742 ierr = VecDestroy(&glee->W);CHKERRQ(ierr); 743 ierr = PetscFree(glee->work);CHKERRQ(ierr); 744 PetscFunctionReturn(0); 745 } 746 747 #undef __FUNCT__ 748 #define __FUNCT__ "TSDestroy_GLEE" 749 static PetscErrorCode TSDestroy_GLEE(TS ts) 750 { 751 PetscErrorCode ierr; 752 753 PetscFunctionBegin; 754 ierr = TSReset_GLEE(ts);CHKERRQ(ierr); 755 ierr = PetscFree(ts->data);CHKERRQ(ierr); 756 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSGLEEGetType_C",NULL);CHKERRQ(ierr); 757 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSGLEESetType_C",NULL);CHKERRQ(ierr); 758 PetscFunctionReturn(0); 759 } 760 761 #undef __FUNCT__ 762 #define __FUNCT__ "TSGLEEGetVecs" 763 static PetscErrorCode TSGLEEGetVecs(TS ts,DM dm,Vec *Ydot) 764 { 765 TS_GLEE *glee = (TS_GLEE*)ts->data; 766 PetscErrorCode ierr; 767 768 PetscFunctionBegin; 769 if (Ydot) { 770 if (dm && dm != ts->dm) { 771 ierr = DMGetNamedGlobalVector(dm,"TSGLEE_Ydot",Ydot);CHKERRQ(ierr); 772 } else *Ydot = glee->Ydot; 773 } 774 PetscFunctionReturn(0); 775 } 776 777 778 #undef __FUNCT__ 779 #define __FUNCT__ "TSGLEERestoreVecs" 780 static PetscErrorCode TSGLEERestoreVecs(TS ts,DM dm,Vec *Ydot) 781 { 782 PetscErrorCode ierr; 783 784 PetscFunctionBegin; 785 if (Ydot) { 786 if (dm && dm != ts->dm) { 787 ierr = DMRestoreNamedGlobalVector(dm,"TSGLEE_Ydot",Ydot);CHKERRQ(ierr); 788 } 789 } 790 PetscFunctionReturn(0); 791 } 792 793 /* 794 This defines the nonlinear equation that is to be solved with SNES 795 */ 796 #undef __FUNCT__ 797 #define __FUNCT__ "SNESTSFormFunction_GLEE" 798 static PetscErrorCode SNESTSFormFunction_GLEE(SNES snes,Vec X,Vec F,TS ts) 799 { 800 TS_GLEE *glee = (TS_GLEE*)ts->data; 801 DM dm,dmsave; 802 Vec Ydot; 803 PetscReal shift = glee->scoeff / ts->time_step; 804 PetscErrorCode ierr; 805 806 PetscFunctionBegin; 807 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 808 ierr = TSGLEEGetVecs(ts,dm,&Ydot);CHKERRQ(ierr); 809 /* Set Ydot = shift*X */ 810 ierr = VecCopy(X,Ydot);CHKERRQ(ierr); 811 ierr = VecScale(Ydot,shift);CHKERRQ(ierr); 812 dmsave = ts->dm; 813 ts->dm = dm; 814 815 ierr = TSComputeIFunction(ts,glee->stage_time,X,Ydot,F,PETSC_FALSE);CHKERRQ(ierr); 816 817 ts->dm = dmsave; 818 ierr = TSGLEERestoreVecs(ts,dm,&Ydot);CHKERRQ(ierr); 819 PetscFunctionReturn(0); 820 } 821 822 #undef __FUNCT__ 823 #define __FUNCT__ "SNESTSFormJacobian_GLEE" 824 static PetscErrorCode SNESTSFormJacobian_GLEE(SNES snes,Vec X,Mat A,Mat B,TS ts) 825 { 826 TS_GLEE *glee = (TS_GLEE*)ts->data; 827 DM dm,dmsave; 828 Vec Ydot; 829 PetscReal shift = glee->scoeff / ts->time_step; 830 PetscErrorCode ierr; 831 832 PetscFunctionBegin; 833 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 834 ierr = TSGLEEGetVecs(ts,dm,&Ydot);CHKERRQ(ierr); 835 /* glee->Ydot has already been computed in SNESTSFormFunction_GLEE (SNES guarantees this) */ 836 dmsave = ts->dm; 837 ts->dm = dm; 838 839 ierr = TSComputeIJacobian(ts,glee->stage_time,X,Ydot,shift,A,B,PETSC_FALSE);CHKERRQ(ierr); 840 841 ts->dm = dmsave; 842 ierr = TSGLEERestoreVecs(ts,dm,&Ydot);CHKERRQ(ierr); 843 PetscFunctionReturn(0); 844 } 845 846 #undef __FUNCT__ 847 #define __FUNCT__ "DMCoarsenHook_TSGLEE" 848 static PetscErrorCode DMCoarsenHook_TSGLEE(DM fine,DM coarse,void *ctx) 849 { 850 PetscFunctionBegin; 851 PetscFunctionReturn(0); 852 } 853 854 #undef __FUNCT__ 855 #define __FUNCT__ "DMRestrictHook_TSGLEE" 856 static PetscErrorCode DMRestrictHook_TSGLEE(DM fine,Mat restrct,Vec rscale,Mat inject,DM coarse,void *ctx) 857 { 858 PetscFunctionBegin; 859 PetscFunctionReturn(0); 860 } 861 862 863 #undef __FUNCT__ 864 #define __FUNCT__ "DMSubDomainHook_TSGLEE" 865 static PetscErrorCode DMSubDomainHook_TSGLEE(DM dm,DM subdm,void *ctx) 866 { 867 PetscFunctionBegin; 868 PetscFunctionReturn(0); 869 } 870 871 #undef __FUNCT__ 872 #define __FUNCT__ "DMSubDomainRestrictHook_TSGLEE" 873 static PetscErrorCode DMSubDomainRestrictHook_TSGLEE(DM dm,VecScatter gscat,VecScatter lscat,DM subdm,void *ctx) 874 { 875 PetscFunctionBegin; 876 PetscFunctionReturn(0); 877 } 878 879 #undef __FUNCT__ 880 #define __FUNCT__ "TSSetUp_GLEE" 881 static PetscErrorCode TSSetUp_GLEE(TS ts) 882 { 883 TS_GLEE *glee = (TS_GLEE*)ts->data; 884 GLEETableau tab; 885 PetscInt s,r; 886 PetscErrorCode ierr; 887 DM dm; 888 889 PetscFunctionBegin; 890 if (!glee->tableau) { 891 ierr = TSGLEESetType(ts,TSGLEEDefaultType);CHKERRQ(ierr); 892 } 893 tab = glee->tableau; 894 s = tab->s; 895 r = tab->r; 896 ierr = VecDuplicateVecs(ts->vec_sol,r,&glee->Y);CHKERRQ(ierr); 897 ierr = VecDuplicateVecs(ts->vec_sol,r,&glee->X);CHKERRQ(ierr); 898 ierr = VecDuplicateVecs(ts->vec_sol,s,&glee->YStage);CHKERRQ(ierr); 899 ierr = VecDuplicateVecs(ts->vec_sol,s,&glee->YdotStage);CHKERRQ(ierr); 900 ierr = VecDuplicate(ts->vec_sol,&glee->Ydot);CHKERRQ(ierr); 901 ierr = VecDuplicate(ts->vec_sol,&glee->yGErr);CHKERRQ(ierr); 902 ierr = VecZeroEntries(glee->yGErr);CHKERRQ(ierr); 903 ierr = VecDuplicate(ts->vec_sol,&glee->W);CHKERRQ(ierr); 904 ierr = PetscMalloc1(s,&glee->work);CHKERRQ(ierr); 905 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 906 if (dm) { 907 ierr = DMCoarsenHookAdd(dm,DMCoarsenHook_TSGLEE,DMRestrictHook_TSGLEE,ts);CHKERRQ(ierr); 908 ierr = DMSubDomainHookAdd(dm,DMSubDomainHook_TSGLEE,DMSubDomainRestrictHook_TSGLEE,ts);CHKERRQ(ierr); 909 } 910 PetscFunctionReturn(0); 911 } 912 913 #undef __FUNCT__ 914 #define __FUNCT__ "TSStartingMethod_GLEE" 915 PetscErrorCode TSStartingMethod_GLEE(TS ts) 916 { 917 TS_GLEE *glee = (TS_GLEE*)ts->data; 918 GLEETableau tab = glee->tableau; 919 PetscInt r=tab->r,i; 920 PetscReal *S=tab->S; 921 PetscErrorCode ierr; 922 923 PetscFunctionBegin; 924 for (i=0; i<r; i++) { 925 ierr = VecZeroEntries(glee->Y[i]);CHKERRQ(ierr); 926 ierr = VecAXPY(glee->Y[i],S[i],ts->vec_sol);CHKERRQ(ierr); 927 } 928 929 PetscFunctionReturn(0); 930 } 931 932 933 /*------------------------------------------------------------*/ 934 935 #undef __FUNCT__ 936 #define __FUNCT__ "TSSetFromOptions_GLEE" 937 static PetscErrorCode TSSetFromOptions_GLEE(PetscOptionItems *PetscOptionsObject,TS ts) 938 { 939 PetscErrorCode ierr; 940 char gleetype[256]; 941 942 PetscFunctionBegin; 943 ierr = PetscOptionsHead(PetscOptionsObject,"GLEE ODE solver options");CHKERRQ(ierr); 944 { 945 GLEETableauLink link; 946 PetscInt count,choice; 947 PetscBool flg; 948 const char **namelist; 949 950 ierr = PetscStrncpy(gleetype,TSGLEEDefaultType,sizeof(gleetype));CHKERRQ(ierr); 951 for (link=GLEETableauList,count=0; link; link=link->next,count++) ; 952 ierr = PetscMalloc1(count,&namelist);CHKERRQ(ierr); 953 for (link=GLEETableauList,count=0; link; link=link->next,count++) namelist[count] = link->tab.name; 954 ierr = PetscOptionsEList("-ts_glee_type","Family of GLEE method","TSGLEESetType",(const char*const*)namelist,count,gleetype,&choice,&flg);CHKERRQ(ierr); 955 ierr = TSGLEESetType(ts,flg ? namelist[choice] : gleetype);CHKERRQ(ierr); 956 ierr = PetscFree(namelist);CHKERRQ(ierr); 957 } 958 ierr = PetscOptionsTail();CHKERRQ(ierr); 959 PetscFunctionReturn(0); 960 } 961 962 #undef __FUNCT__ 963 #define __FUNCT__ "PetscFormatRealArray" 964 static PetscErrorCode PetscFormatRealArray(char buf[],size_t len,const char *fmt,PetscInt n,const PetscReal x[]) 965 { 966 PetscErrorCode ierr; 967 PetscInt i; 968 size_t left,count; 969 char *p; 970 971 PetscFunctionBegin; 972 for (i=0,p=buf,left=len; i<n; i++) { 973 ierr = PetscSNPrintfCount(p,left,fmt,&count,x[i]);CHKERRQ(ierr); 974 if (count >= left) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Insufficient space in buffer"); 975 left -= count; 976 p += count; 977 *p++ = ' '; 978 } 979 p[i ? 0 : -1] = 0; 980 PetscFunctionReturn(0); 981 } 982 983 #undef __FUNCT__ 984 #define __FUNCT__ "TSView_GLEE" 985 static PetscErrorCode TSView_GLEE(TS ts,PetscViewer viewer) 986 { 987 TS_GLEE *glee = (TS_GLEE*)ts->data; 988 GLEETableau tab = glee->tableau; 989 PetscBool iascii; 990 PetscErrorCode ierr; 991 TSAdapt adapt; 992 993 PetscFunctionBegin; 994 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 995 if (iascii) { 996 TSGLEEType gleetype; 997 char buf[512]; 998 ierr = TSGLEEGetType(ts,&gleetype);CHKERRQ(ierr); 999 ierr = PetscViewerASCIIPrintf(viewer," GLEE %s\n",gleetype);CHKERRQ(ierr); 1000 ierr = PetscFormatRealArray(buf,sizeof(buf),"% 8.6f",tab->s,tab->c);CHKERRQ(ierr); 1001 ierr = PetscViewerASCIIPrintf(viewer," Abscissa c = %s\n",buf);CHKERRQ(ierr); 1002 /* Note: print out r as well */ 1003 } 1004 ierr = TSGetAdapt(ts,&adapt);CHKERRQ(ierr); 1005 ierr = TSAdaptView(adapt,viewer);CHKERRQ(ierr); 1006 PetscFunctionReturn(0); 1007 } 1008 1009 #undef __FUNCT__ 1010 #define __FUNCT__ "TSLoad_GLEE" 1011 static PetscErrorCode TSLoad_GLEE(TS ts,PetscViewer viewer) 1012 { 1013 PetscErrorCode ierr; 1014 SNES snes; 1015 TSAdapt tsadapt; 1016 1017 PetscFunctionBegin; 1018 ierr = TSGetAdapt(ts,&tsadapt);CHKERRQ(ierr); 1019 ierr = TSAdaptLoad(tsadapt,viewer);CHKERRQ(ierr); 1020 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 1021 ierr = SNESLoad(snes,viewer);CHKERRQ(ierr); 1022 /* function and Jacobian context for SNES when used with TS is always ts object */ 1023 ierr = SNESSetFunction(snes,NULL,NULL,ts);CHKERRQ(ierr); 1024 ierr = SNESSetJacobian(snes,NULL,NULL,NULL,ts);CHKERRQ(ierr); 1025 PetscFunctionReturn(0); 1026 } 1027 1028 #undef __FUNCT__ 1029 #define __FUNCT__ "TSGLEESetType" 1030 /*@C 1031 TSGLEESetType - Set the type of GLEE scheme 1032 1033 Logically collective 1034 1035 Input Parameter: 1036 + ts - timestepping context 1037 - gleetype - type of GLEE-scheme 1038 1039 Level: intermediate 1040 1041 .seealso: TSGLEEGetType(), TSGLEE 1042 @*/ 1043 PetscErrorCode TSGLEESetType(TS ts,TSGLEEType gleetype) 1044 { 1045 PetscErrorCode ierr; 1046 1047 PetscFunctionBegin; 1048 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1049 ierr = PetscTryMethod(ts,"TSGLEESetType_C",(TS,TSGLEEType),(ts,gleetype));CHKERRQ(ierr); 1050 PetscFunctionReturn(0); 1051 } 1052 1053 #undef __FUNCT__ 1054 #define __FUNCT__ "TSGLEEGetType" 1055 /*@C 1056 TSGLEEGetType - Get the type of GLEE scheme 1057 1058 Logically collective 1059 1060 Input Parameter: 1061 . ts - timestepping context 1062 1063 Output Parameter: 1064 . gleetype - type of GLEE-scheme 1065 1066 Level: intermediate 1067 1068 .seealso: TSGLEESetType() 1069 @*/ 1070 PetscErrorCode TSGLEEGetType(TS ts,TSGLEEType *gleetype) 1071 { 1072 PetscErrorCode ierr; 1073 1074 PetscFunctionBegin; 1075 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1076 ierr = PetscUseMethod(ts,"TSGLEEGetType_C",(TS,TSGLEEType*),(ts,gleetype));CHKERRQ(ierr); 1077 PetscFunctionReturn(0); 1078 } 1079 1080 #undef __FUNCT__ 1081 #define __FUNCT__ "TSGLEEGetType_GLEE" 1082 PetscErrorCode TSGLEEGetType_GLEE(TS ts,TSGLEEType *gleetype) 1083 { 1084 TS_GLEE *glee = (TS_GLEE*)ts->data; 1085 PetscErrorCode ierr; 1086 1087 PetscFunctionBegin; 1088 if (!glee->tableau) { 1089 ierr = TSGLEESetType(ts,TSGLEEDefaultType);CHKERRQ(ierr); 1090 } 1091 *gleetype = glee->tableau->name; 1092 PetscFunctionReturn(0); 1093 } 1094 #undef __FUNCT__ 1095 #define __FUNCT__ "TSGLEESetType_GLEE" 1096 PetscErrorCode TSGLEESetType_GLEE(TS ts,TSGLEEType gleetype) 1097 { 1098 TS_GLEE *glee = (TS_GLEE*)ts->data; 1099 PetscErrorCode ierr; 1100 PetscBool match; 1101 GLEETableauLink link; 1102 1103 PetscFunctionBegin; 1104 if (glee->tableau) { 1105 ierr = PetscStrcmp(glee->tableau->name,gleetype,&match);CHKERRQ(ierr); 1106 if (match) PetscFunctionReturn(0); 1107 } 1108 for (link = GLEETableauList; link; link=link->next) { 1109 ierr = PetscStrcmp(link->tab.name,gleetype,&match);CHKERRQ(ierr); 1110 if (match) { 1111 ierr = TSReset_GLEE(ts);CHKERRQ(ierr); 1112 glee->tableau = &link->tab; 1113 PetscFunctionReturn(0); 1114 } 1115 } 1116 SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_UNKNOWN_TYPE,"Could not find '%s'",gleetype); 1117 PetscFunctionReturn(0); 1118 } 1119 1120 #undef __FUNCT__ 1121 #define __FUNCT__ "TSGetStages_GLEE" 1122 static PetscErrorCode TSGetStages_GLEE(TS ts,PetscInt *ns,Vec **Y) 1123 { 1124 TS_GLEE *glee = (TS_GLEE*)ts->data; 1125 1126 PetscFunctionBegin; 1127 *ns = glee->tableau->s; 1128 if(Y) *Y = glee->Y; 1129 PetscFunctionReturn(0); 1130 } 1131 1132 #undef __FUNCT__ 1133 #define __FUNCT__ "TSGetSolutionComponents_GLEE" 1134 PetscErrorCode TSGetSolutionComponents_GLEE(TS ts,PetscInt *n,Vec *Y) 1135 { 1136 TS_GLEE *glee = (TS_GLEE*)ts->data; 1137 GLEETableau tab = glee->tableau; 1138 PetscErrorCode ierr; 1139 1140 PetscFunctionBegin; 1141 if (!Y) *n = tab->r; 1142 else { 1143 if ((*n >= 0) && (*n < tab->r)) { 1144 ierr = VecCopy(glee->Y[*n],*Y); CHKERRQ(ierr); 1145 } else SETERRQ3(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_OUTOFRANGE,"Second argument (%d) out of range[%d,%d].",*n,0,tab->r-1); 1146 } 1147 PetscFunctionReturn(0); 1148 } 1149 1150 #undef __FUNCT__ 1151 #define __FUNCT__ "TSGetAuxSolution_GLEE" 1152 PetscErrorCode TSGetAuxSolution_GLEE(TS ts,Vec *X) 1153 { 1154 TS_GLEE *glee = (TS_GLEE*)ts->data; 1155 GLEETableau tab = glee->tableau; 1156 PetscReal *F = tab->Fembed; 1157 PetscInt r = tab->r; 1158 Vec *Y = glee->Y; 1159 PetscErrorCode ierr; 1160 1161 PetscFunctionBegin; 1162 ierr = VecZeroEntries(*X);CHKERRQ(ierr); 1163 ierr = VecMAXPY((*X),r,F,Y);CHKERRQ(ierr); 1164 PetscFunctionReturn(0); 1165 } 1166 1167 #undef __FUNCT__ 1168 #define __FUNCT__ "TSGetTimeError_GLEE" 1169 PetscErrorCode TSGetTimeError_GLEE(TS ts,PetscInt n,Vec *X) 1170 { 1171 TS_GLEE *glee = (TS_GLEE*)ts->data; 1172 GLEETableau tab = glee->tableau; 1173 PetscReal *F = tab->Ferror; 1174 PetscInt r = tab->r; 1175 Vec *Y = glee->Y; 1176 PetscErrorCode ierr; 1177 1178 PetscFunctionBegin; 1179 ierr = VecZeroEntries(*X);CHKERRQ(ierr); 1180 if(n==0){ 1181 ierr = VecMAXPY((*X),r,F,Y);CHKERRQ(ierr); 1182 } else if(n==-1) { 1183 *X=glee->yGErr; 1184 } 1185 PetscFunctionReturn(0); 1186 } 1187 1188 #undef __FUNCT__ 1189 #define __FUNCT__ "TSSetTimeError_GLEE" 1190 PetscErrorCode TSSetTimeError_GLEE(TS ts,Vec X) 1191 { 1192 TS_GLEE *glee = (TS_GLEE*)ts->data; 1193 GLEETableau tab = glee->tableau; 1194 PetscReal *S = tab->Serror; 1195 PetscInt r = tab->r,i; 1196 Vec *Y = glee->Y; 1197 PetscErrorCode ierr; 1198 1199 PetscFunctionBegin; 1200 if (r != 2) SETERRQ2(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"TSSetTimeError_GLEE not supported for '%s' with r=%D.",tab->name,tab->r); 1201 for (i=1; i<r; i++) { 1202 ierr = VecCopy(ts->vec_sol,Y[i]); CHKERRQ(ierr); 1203 ierr = VecAXPBY(Y[i],S[0],S[1],X); CHKERRQ(ierr); 1204 ierr = VecCopy(X,glee->yGErr); CHKERRQ(ierr); 1205 } 1206 PetscFunctionReturn(0); 1207 } 1208 1209 /* ------------------------------------------------------------ */ 1210 /*MC 1211 TSGLEE - ODE and DAE solver using General Linear with Error Estimation schemes 1212 1213 The user should provide the right hand side of the equation 1214 using TSSetRHSFunction(). 1215 1216 Notes: 1217 The default is TSGLEE35, it can be changed with TSGLEESetType() or -ts_glee_type 1218 1219 Level: beginner 1220 1221 .seealso: TSCreate(), TS, TSSetType(), TSGLEESetType(), TSGLEEGetType(), 1222 TSGLEE23, TTSGLEE24, TSGLEE35, TSGLEE25I, TSGLEEEXRK2A, 1223 TSGLEERK32G1, TSGLEERK285EX, TSGLEEType, TSGLEERegister() 1224 1225 M*/ 1226 #undef __FUNCT__ 1227 #define __FUNCT__ "TSCreate_GLEE" 1228 PETSC_EXTERN PetscErrorCode TSCreate_GLEE(TS ts) 1229 { 1230 TS_GLEE *th; 1231 PetscErrorCode ierr; 1232 1233 PetscFunctionBegin; 1234 #if !defined(PETSC_USE_DYNAMIC_LIBRARIES) 1235 ierr = TSGLEEInitializePackage();CHKERRQ(ierr); 1236 #endif 1237 1238 ts->ops->reset = TSReset_GLEE; 1239 ts->ops->destroy = TSDestroy_GLEE; 1240 ts->ops->view = TSView_GLEE; 1241 ts->ops->load = TSLoad_GLEE; 1242 ts->ops->setup = TSSetUp_GLEE; 1243 ts->ops->step = TSStep_GLEE; 1244 ts->ops->interpolate = TSInterpolate_GLEE; 1245 ts->ops->evaluatestep = TSEvaluateStep_GLEE; 1246 ts->ops->setfromoptions = TSSetFromOptions_GLEE; 1247 ts->ops->getstages = TSGetStages_GLEE; 1248 ts->ops->snesfunction = SNESTSFormFunction_GLEE; 1249 ts->ops->snesjacobian = SNESTSFormJacobian_GLEE; 1250 ts->ops->getsolutioncomponents = TSGetSolutionComponents_GLEE; 1251 ts->ops->getauxsolution = TSGetAuxSolution_GLEE; 1252 ts->ops->gettimeerror = TSGetTimeError_GLEE; 1253 ts->ops->settimeerror = TSSetTimeError_GLEE; 1254 ts->ops->startingmethod = TSStartingMethod_GLEE; 1255 1256 ierr = PetscNewLog(ts,&th);CHKERRQ(ierr); 1257 ts->data = (void*)th; 1258 1259 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSGLEEGetType_C",TSGLEEGetType_GLEE);CHKERRQ(ierr); 1260 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSGLEESetType_C",TSGLEESetType_GLEE);CHKERRQ(ierr); 1261 PetscFunctionReturn(0); 1262 } 1263