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