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 TSInitializePackage(). 395 396 Level: developer 397 398 .keywords: TS, TSGLEE, initialize, package 399 .seealso: PetscInitialize() 400 @*/ 401 PetscErrorCode TSGLEEInitializePackage(void) 402 { 403 PetscErrorCode ierr; 404 405 PetscFunctionBegin; 406 if (TSGLEEPackageInitialized) PetscFunctionReturn(0); 407 TSGLEEPackageInitialized = PETSC_TRUE; 408 ierr = TSGLEERegisterAll();CHKERRQ(ierr); 409 ierr = PetscObjectComposedDataRegister(&explicit_stage_time_id);CHKERRQ(ierr); 410 ierr = PetscRegisterFinalize(TSGLEEFinalizePackage);CHKERRQ(ierr); 411 PetscFunctionReturn(0); 412 } 413 414 /*@C 415 TSGLEEFinalizePackage - This function destroys everything in the TSGLEE package. It is 416 called from PetscFinalize(). 417 418 Level: developer 419 420 .keywords: Petsc, destroy, package 421 .seealso: PetscFinalize() 422 @*/ 423 PetscErrorCode TSGLEEFinalizePackage(void) 424 { 425 PetscErrorCode ierr; 426 427 PetscFunctionBegin; 428 TSGLEEPackageInitialized = PETSC_FALSE; 429 ierr = TSGLEERegisterDestroy();CHKERRQ(ierr); 430 PetscFunctionReturn(0); 431 } 432 433 /*@C 434 TSGLEERegister - register an GLEE scheme by providing the entries in the Butcher tableau 435 436 Not Collective, but the same schemes should be registered on all processes on which they will be used 437 438 Input Parameters: 439 + name - identifier for method 440 . order - order of method 441 . s - number of stages 442 . r - number of steps 443 . gamma - LTE ratio 444 . A - stage coefficients (dimension s*s, row-major) 445 . B - step completion coefficients (dimension r*s, row-major) 446 . U - method coefficients (dimension s*r, row-major) 447 . V - method coefficients (dimension r*r, row-major) 448 . S - starting coefficients 449 . F - finishing coefficients 450 . c - abscissa (dimension s; NULL to use row sums of A) 451 . Fembed - step completion coefficients for embedded method 452 . Ferror - error computation coefficients 453 . Serror - error initialization coefficients 454 . pinterp - order of interpolation (0 if unavailable) 455 - binterp - array of interpolation coefficients (NULL if unavailable) 456 457 Notes: 458 Several GLEE methods are provided, this function is only needed to create new methods. 459 460 Level: advanced 461 462 .keywords: TS, register 463 464 .seealso: TSGLEE 465 @*/ 466 PetscErrorCode TSGLEERegister(TSGLEEType name,PetscInt order,PetscInt s, PetscInt r, 467 PetscReal gamma, 468 const PetscReal A[],const PetscReal B[], 469 const PetscReal U[],const PetscReal V[], 470 const PetscReal S[],const PetscReal F[], 471 const PetscReal c[], 472 const PetscReal Fembed[],const PetscReal Ferror[], 473 const PetscReal Serror[], 474 PetscInt pinterp, const PetscReal binterp[]) 475 { 476 PetscErrorCode ierr; 477 GLEETableauLink link; 478 GLEETableau t; 479 PetscInt i,j; 480 481 PetscFunctionBegin; 482 ierr = TSGLEEInitializePackage();CHKERRQ(ierr); 483 ierr = PetscMalloc(sizeof(*link),&link);CHKERRQ(ierr); 484 ierr = PetscMemzero(link,sizeof(*link));CHKERRQ(ierr); 485 t = &link->tab; 486 ierr = PetscStrallocpy(name,&t->name);CHKERRQ(ierr); 487 t->order = order; 488 t->s = s; 489 t->r = r; 490 t->gamma = gamma; 491 ierr = PetscMalloc5(s*s,&t->A,r*r,&t->V,s,&t->c,r*s,&t->B,s*r,&t->U); CHKERRQ(ierr); 492 ierr = PetscMalloc2(r,&t->S,r,&t->F); CHKERRQ(ierr); 493 ierr = PetscMemcpy(t->A,A,s*s*sizeof(A[0]));CHKERRQ(ierr); 494 ierr = PetscMemcpy(t->B,B,r*s*sizeof(B[0]));CHKERRQ(ierr); 495 ierr = PetscMemcpy(t->U,U,s*r*sizeof(U[0]));CHKERRQ(ierr); 496 ierr = PetscMemcpy(t->V,V,r*r*sizeof(V[0]));CHKERRQ(ierr); 497 ierr = PetscMemcpy(t->S,S,r *sizeof(S[0]));CHKERRQ(ierr); 498 ierr = PetscMemcpy(t->F,F,r *sizeof(F[0]));CHKERRQ(ierr); 499 if (c) { 500 ierr = PetscMemcpy(t->c,c,s*sizeof(c[0])); CHKERRQ(ierr); 501 } else { 502 for (i=0; i<s; i++) for (j=0,t->c[i]=0; j<s; j++) t->c[i] += A[i*s+j]; 503 } 504 ierr = PetscMalloc1(r,&t->Fembed);CHKERRQ(ierr); 505 ierr = PetscMalloc1(r,&t->Ferror);CHKERRQ(ierr); 506 ierr = PetscMalloc1(r,&t->Serror);CHKERRQ(ierr); 507 ierr = PetscMemcpy(t->Fembed,Fembed,r*sizeof(Fembed[0]));CHKERRQ(ierr); 508 ierr = PetscMemcpy(t->Ferror,Ferror,r*sizeof(Ferror[0]));CHKERRQ(ierr); 509 ierr = PetscMemcpy(t->Serror,Serror,r*sizeof(Serror[0]));CHKERRQ(ierr); 510 t->pinterp = pinterp; 511 ierr = PetscMalloc1(s*pinterp,&t->binterp);CHKERRQ(ierr); 512 ierr = PetscMemcpy(t->binterp,binterp,s*pinterp*sizeof(binterp[0]));CHKERRQ(ierr); 513 514 link->next = GLEETableauList; 515 GLEETableauList = link; 516 PetscFunctionReturn(0); 517 } 518 519 static PetscErrorCode TSEvaluateStep_GLEE(TS ts,PetscInt order,Vec X,PetscBool *done) 520 { 521 TS_GLEE *glee = (TS_GLEE*) ts->data; 522 GLEETableau tab = glee->tableau; 523 PetscReal h, *B = tab->B, *V = tab->V, 524 *F = tab->F, 525 *Fembed = tab->Fembed; 526 PetscInt s = tab->s, r = tab->r, i, j; 527 Vec *Y = glee->Y, *YdotStage = glee->YdotStage; 528 PetscScalar *ws = glee->swork, *wr = glee->rwork; 529 PetscErrorCode ierr; 530 531 PetscFunctionBegin; 532 533 switch (glee->status) { 534 case TS_STEP_INCOMPLETE: 535 case TS_STEP_PENDING: 536 h = ts->time_step; break; 537 case TS_STEP_COMPLETE: 538 h = ts->ptime - ts->ptime_prev; break; 539 default: SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_PLIB,"Invalid TSStepStatus"); 540 } 541 542 if (order == tab->order) { 543 544 /* Note: Irrespective of whether status is TS_STEP_INCOMPLETE 545 or TS_STEP_COMPLETE, glee->X has the solution at the 546 beginning of the time step. So no need to roll-back. 547 */ 548 if (glee->status == TS_STEP_INCOMPLETE) { 549 for (i=0; i<r; i++) { 550 ierr = VecZeroEntries(Y[i]); CHKERRQ(ierr); 551 for (j=0; j<r; j++) wr[j] = V[i*r+j]; 552 ierr = VecMAXPY(Y[i],r,wr,glee->X); CHKERRQ(ierr); 553 for (j=0; j<s; j++) ws[j] = h*B[i*s+j]; 554 ierr = VecMAXPY(Y[i],s,ws,YdotStage); CHKERRQ(ierr); 555 } 556 ierr = VecZeroEntries(X);CHKERRQ(ierr); 557 for (j=0; j<r; j++) wr[j] = F[j]; 558 ierr = VecMAXPY(X,r,wr,Y);CHKERRQ(ierr); 559 } else {ierr = VecCopy(ts->vec_sol,X);CHKERRQ(ierr);} 560 PetscFunctionReturn(0); 561 562 } else if (order == tab->order-1) { 563 564 /* Complete with the embedded method (Fembed) */ 565 for (i=0; i<r; i++) { 566 ierr = VecZeroEntries(Y[i]); CHKERRQ(ierr); 567 for (j=0; j<r; j++) wr[j] = V[i*r+j]; 568 ierr = VecMAXPY(Y[i],r,wr,glee->X); CHKERRQ(ierr); 569 for (j=0; j<s; j++) ws[j] = h*B[i*s+j]; 570 ierr = VecMAXPY(Y[i],s,ws,YdotStage); CHKERRQ(ierr); 571 } 572 ierr = VecZeroEntries(X);CHKERRQ(ierr); 573 for (j=0; j<r; j++) wr[j] = Fembed[j]; 574 ierr = VecMAXPY(X,r,wr,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 static PetscErrorCode TSStep_GLEE(TS ts) 585 { 586 TS_GLEE *glee = (TS_GLEE*)ts->data; 587 GLEETableau tab = glee->tableau; 588 const PetscInt s = tab->s, r = tab->r; 589 PetscReal *A = tab->A, *U = tab->U, 590 *F = tab->F, 591 *c = tab->c; 592 Vec *Y = glee->Y, *X = glee->X, 593 *YStage = glee->YStage, 594 *YdotStage = glee->YdotStage, 595 W = glee->W; 596 SNES snes; 597 PetscScalar *ws = glee->swork, *wr = glee->rwork; 598 TSAdapt adapt; 599 PetscInt i,j,reject,next_scheme,its,lits; 600 PetscReal next_time_step; 601 PetscReal t; 602 PetscBool accept; 603 PetscErrorCode ierr; 604 605 PetscFunctionBegin; 606 ierr = PetscCitationsRegister(citation,&cited);CHKERRQ(ierr); 607 608 for (i=0; i<r; i++) { ierr = VecCopy(Y[i],X[i]); CHKERRQ(ierr); } 609 610 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 611 next_time_step = ts->time_step; 612 t = ts->ptime; 613 accept = PETSC_TRUE; 614 glee->status = TS_STEP_INCOMPLETE; 615 616 for (reject=0; reject<ts->max_reject && !ts->reason; reject++,ts->reject++) { 617 618 PetscReal h = ts->time_step; 619 ierr = TSPreStep(ts);CHKERRQ(ierr); 620 621 for (i=0; i<s; i++) { 622 623 glee->stage_time = t + h*c[i]; 624 ierr = TSPreStage(ts,glee->stage_time); CHKERRQ(ierr); 625 626 if (A[i*s+i] == 0) { /* Explicit stage */ 627 ierr = VecZeroEntries(YStage[i]);CHKERRQ(ierr); 628 for (j=0; j<r; j++) wr[j] = U[i*r+j]; 629 ierr = VecMAXPY(YStage[i],r,wr,X);CHKERRQ(ierr); 630 for (j=0; j<i; j++) ws[j] = h*A[i*s+j]; 631 ierr = VecMAXPY(YStage[i],i,ws,YdotStage);CHKERRQ(ierr); 632 } else { /* Implicit stage */ 633 glee->scoeff = 1.0/A[i*s+i]; 634 /* compute right-hand-side */ 635 ierr = VecZeroEntries(W);CHKERRQ(ierr); 636 for (j=0; j<r; j++) wr[j] = U[i*r+j]; 637 ierr = VecMAXPY(W,r,wr,X);CHKERRQ(ierr); 638 for (j=0; j<i; j++) ws[j] = h*A[i*s+j]; 639 ierr = VecMAXPY(W,i,ws,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,(PetscReal)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 for (j=0; j<r; j++) wr[j] = F[j]; 675 ierr = VecMAXPY(ts->vec_sol,r,wr,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 static PetscErrorCode TSInterpolate_GLEE(TS ts,PetscReal itime,Vec X) 686 { 687 TS_GLEE *glee = (TS_GLEE*)ts->data; 688 PetscInt s=glee->tableau->s, pinterp=glee->tableau->pinterp,i,j; 689 PetscReal h,tt,t; 690 PetscScalar *b; 691 const PetscReal *B = glee->tableau->binterp; 692 PetscErrorCode ierr; 693 694 PetscFunctionBegin; 695 if (!B) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"TSGLEE %s does not have an interpolation formula",glee->tableau->name); 696 switch (glee->status) { 697 case TS_STEP_INCOMPLETE: 698 case TS_STEP_PENDING: 699 h = ts->time_step; 700 t = (itime - ts->ptime)/h; 701 break; 702 case TS_STEP_COMPLETE: 703 h = ts->ptime - ts->ptime_prev; 704 t = (itime - ts->ptime)/h + 1; /* In the interval [0,1] */ 705 break; 706 default: SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_PLIB,"Invalid TSStepStatus"); 707 } 708 ierr = PetscMalloc1(s,&b);CHKERRQ(ierr); 709 for (i=0; i<s; i++) b[i] = 0; 710 for (j=0,tt=t; j<pinterp; j++,tt*=t) { 711 for (i=0; i<s; i++) { 712 b[i] += h * B[i*pinterp+j] * tt; 713 } 714 } 715 ierr = VecCopy(glee->YStage[0],X);CHKERRQ(ierr); 716 ierr = VecMAXPY(X,s,b,glee->YdotStage);CHKERRQ(ierr); 717 ierr = PetscFree(b);CHKERRQ(ierr); 718 PetscFunctionReturn(0); 719 } 720 721 /*------------------------------------------------------------*/ 722 static PetscErrorCode TSReset_GLEE(TS ts) 723 { 724 TS_GLEE *glee = (TS_GLEE*)ts->data; 725 PetscInt s, r; 726 PetscErrorCode ierr; 727 728 PetscFunctionBegin; 729 if (!glee->tableau) PetscFunctionReturn(0); 730 s = glee->tableau->s; 731 r = glee->tableau->r; 732 ierr = VecDestroyVecs(r,&glee->Y);CHKERRQ(ierr); 733 ierr = VecDestroyVecs(r,&glee->X);CHKERRQ(ierr); 734 ierr = VecDestroyVecs(s,&glee->YStage);CHKERRQ(ierr); 735 ierr = VecDestroyVecs(s,&glee->YdotStage);CHKERRQ(ierr); 736 ierr = VecDestroy(&glee->Ydot);CHKERRQ(ierr); 737 ierr = VecDestroy(&glee->yGErr);CHKERRQ(ierr); 738 ierr = VecDestroy(&glee->W);CHKERRQ(ierr); 739 ierr = PetscFree2(glee->swork,glee->rwork);CHKERRQ(ierr); 740 PetscFunctionReturn(0); 741 } 742 743 static PetscErrorCode TSGLEEGetVecs(TS ts,DM dm,Vec *Ydot) 744 { 745 TS_GLEE *glee = (TS_GLEE*)ts->data; 746 PetscErrorCode ierr; 747 748 PetscFunctionBegin; 749 if (Ydot) { 750 if (dm && dm != ts->dm) { 751 ierr = DMGetNamedGlobalVector(dm,"TSGLEE_Ydot",Ydot);CHKERRQ(ierr); 752 } else *Ydot = glee->Ydot; 753 } 754 PetscFunctionReturn(0); 755 } 756 757 758 static PetscErrorCode TSGLEERestoreVecs(TS ts,DM dm,Vec *Ydot) 759 { 760 PetscErrorCode ierr; 761 762 PetscFunctionBegin; 763 if (Ydot) { 764 if (dm && dm != ts->dm) { 765 ierr = DMRestoreNamedGlobalVector(dm,"TSGLEE_Ydot",Ydot);CHKERRQ(ierr); 766 } 767 } 768 PetscFunctionReturn(0); 769 } 770 771 /* 772 This defines the nonlinear equation that is to be solved with SNES 773 */ 774 static PetscErrorCode SNESTSFormFunction_GLEE(SNES snes,Vec X,Vec F,TS ts) 775 { 776 TS_GLEE *glee = (TS_GLEE*)ts->data; 777 DM dm,dmsave; 778 Vec Ydot; 779 PetscReal shift = glee->scoeff / ts->time_step; 780 PetscErrorCode ierr; 781 782 PetscFunctionBegin; 783 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 784 ierr = TSGLEEGetVecs(ts,dm,&Ydot);CHKERRQ(ierr); 785 /* Set Ydot = shift*X */ 786 ierr = VecCopy(X,Ydot);CHKERRQ(ierr); 787 ierr = VecScale(Ydot,shift);CHKERRQ(ierr); 788 dmsave = ts->dm; 789 ts->dm = dm; 790 791 ierr = TSComputeIFunction(ts,glee->stage_time,X,Ydot,F,PETSC_FALSE);CHKERRQ(ierr); 792 793 ts->dm = dmsave; 794 ierr = TSGLEERestoreVecs(ts,dm,&Ydot);CHKERRQ(ierr); 795 PetscFunctionReturn(0); 796 } 797 798 static PetscErrorCode SNESTSFormJacobian_GLEE(SNES snes,Vec X,Mat A,Mat B,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 /* glee->Ydot has already been computed in SNESTSFormFunction_GLEE (SNES guarantees this) */ 810 dmsave = ts->dm; 811 ts->dm = dm; 812 813 ierr = TSComputeIJacobian(ts,glee->stage_time,X,Ydot,shift,A,B,PETSC_FALSE);CHKERRQ(ierr); 814 815 ts->dm = dmsave; 816 ierr = TSGLEERestoreVecs(ts,dm,&Ydot);CHKERRQ(ierr); 817 PetscFunctionReturn(0); 818 } 819 820 static PetscErrorCode DMCoarsenHook_TSGLEE(DM fine,DM coarse,void *ctx) 821 { 822 PetscFunctionBegin; 823 PetscFunctionReturn(0); 824 } 825 826 static PetscErrorCode DMRestrictHook_TSGLEE(DM fine,Mat restrct,Vec rscale,Mat inject,DM coarse,void *ctx) 827 { 828 PetscFunctionBegin; 829 PetscFunctionReturn(0); 830 } 831 832 833 static PetscErrorCode DMSubDomainHook_TSGLEE(DM dm,DM subdm,void *ctx) 834 { 835 PetscFunctionBegin; 836 PetscFunctionReturn(0); 837 } 838 839 static PetscErrorCode DMSubDomainRestrictHook_TSGLEE(DM dm,VecScatter gscat,VecScatter lscat,DM subdm,void *ctx) 840 { 841 PetscFunctionBegin; 842 PetscFunctionReturn(0); 843 } 844 845 static PetscErrorCode TSSetUp_GLEE(TS ts) 846 { 847 TS_GLEE *glee = (TS_GLEE*)ts->data; 848 GLEETableau tab; 849 PetscInt s,r; 850 PetscErrorCode ierr; 851 DM dm; 852 853 PetscFunctionBegin; 854 if (!glee->tableau) { 855 ierr = TSGLEESetType(ts,TSGLEEDefaultType);CHKERRQ(ierr); 856 } 857 tab = glee->tableau; 858 s = tab->s; 859 r = tab->r; 860 ierr = VecDuplicateVecs(ts->vec_sol,r,&glee->Y);CHKERRQ(ierr); 861 ierr = VecDuplicateVecs(ts->vec_sol,r,&glee->X);CHKERRQ(ierr); 862 ierr = VecDuplicateVecs(ts->vec_sol,s,&glee->YStage);CHKERRQ(ierr); 863 ierr = VecDuplicateVecs(ts->vec_sol,s,&glee->YdotStage);CHKERRQ(ierr); 864 ierr = VecDuplicate(ts->vec_sol,&glee->Ydot);CHKERRQ(ierr); 865 ierr = VecDuplicate(ts->vec_sol,&glee->yGErr);CHKERRQ(ierr); 866 ierr = VecZeroEntries(glee->yGErr);CHKERRQ(ierr); 867 ierr = VecDuplicate(ts->vec_sol,&glee->W);CHKERRQ(ierr); 868 ierr = PetscMalloc2(s,&glee->swork,r,&glee->rwork); CHKERRQ(ierr); 869 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 870 ierr = DMCoarsenHookAdd(dm,DMCoarsenHook_TSGLEE,DMRestrictHook_TSGLEE,ts);CHKERRQ(ierr); 871 ierr = DMSubDomainHookAdd(dm,DMSubDomainHook_TSGLEE,DMSubDomainRestrictHook_TSGLEE,ts);CHKERRQ(ierr); 872 PetscFunctionReturn(0); 873 } 874 875 PetscErrorCode TSStartingMethod_GLEE(TS ts) 876 { 877 TS_GLEE *glee = (TS_GLEE*)ts->data; 878 GLEETableau tab = glee->tableau; 879 PetscInt r=tab->r,i; 880 PetscReal *S=tab->S; 881 PetscErrorCode ierr; 882 883 PetscFunctionBegin; 884 for (i=0; i<r; i++) { 885 ierr = VecZeroEntries(glee->Y[i]);CHKERRQ(ierr); 886 ierr = VecAXPY(glee->Y[i],S[i],ts->vec_sol);CHKERRQ(ierr); 887 } 888 889 PetscFunctionReturn(0); 890 } 891 892 893 /*------------------------------------------------------------*/ 894 895 static PetscErrorCode TSSetFromOptions_GLEE(PetscOptionItems *PetscOptionsObject,TS ts) 896 { 897 PetscErrorCode ierr; 898 char gleetype[256]; 899 900 PetscFunctionBegin; 901 ierr = PetscOptionsHead(PetscOptionsObject,"GLEE ODE solver options");CHKERRQ(ierr); 902 { 903 GLEETableauLink link; 904 PetscInt count,choice; 905 PetscBool flg; 906 const char **namelist; 907 908 ierr = PetscStrncpy(gleetype,TSGLEEDefaultType,sizeof(gleetype));CHKERRQ(ierr); 909 for (link=GLEETableauList,count=0; link; link=link->next,count++) ; 910 ierr = PetscMalloc1(count,(char***)&namelist);CHKERRQ(ierr); 911 for (link=GLEETableauList,count=0; link; link=link->next,count++) namelist[count] = link->tab.name; 912 ierr = PetscOptionsEList("-ts_glee_type","Family of GLEE method","TSGLEESetType",(const char*const*)namelist,count,gleetype,&choice,&flg);CHKERRQ(ierr); 913 ierr = TSGLEESetType(ts,flg ? namelist[choice] : gleetype);CHKERRQ(ierr); 914 ierr = PetscFree(namelist);CHKERRQ(ierr); 915 } 916 ierr = PetscOptionsTail();CHKERRQ(ierr); 917 PetscFunctionReturn(0); 918 } 919 920 static PetscErrorCode TSView_GLEE(TS ts,PetscViewer viewer) 921 { 922 TS_GLEE *glee = (TS_GLEE*)ts->data; 923 GLEETableau tab = glee->tableau; 924 PetscBool iascii; 925 PetscErrorCode ierr; 926 927 PetscFunctionBegin; 928 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 929 if (iascii) { 930 TSGLEEType gleetype; 931 char buf[512]; 932 ierr = TSGLEEGetType(ts,&gleetype);CHKERRQ(ierr); 933 ierr = PetscViewerASCIIPrintf(viewer," GLEE type %s\n",gleetype);CHKERRQ(ierr); 934 ierr = PetscFormatRealArray(buf,sizeof(buf),"% 8.6f",tab->s,tab->c);CHKERRQ(ierr); 935 ierr = PetscViewerASCIIPrintf(viewer," Abscissa c = %s\n",buf);CHKERRQ(ierr); 936 /* Note: print out r as well */ 937 } 938 PetscFunctionReturn(0); 939 } 940 941 static PetscErrorCode TSLoad_GLEE(TS ts,PetscViewer viewer) 942 { 943 PetscErrorCode ierr; 944 SNES snes; 945 TSAdapt tsadapt; 946 947 PetscFunctionBegin; 948 ierr = TSGetAdapt(ts,&tsadapt);CHKERRQ(ierr); 949 ierr = TSAdaptLoad(tsadapt,viewer);CHKERRQ(ierr); 950 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 951 ierr = SNESLoad(snes,viewer);CHKERRQ(ierr); 952 /* function and Jacobian context for SNES when used with TS is always ts object */ 953 ierr = SNESSetFunction(snes,NULL,NULL,ts);CHKERRQ(ierr); 954 ierr = SNESSetJacobian(snes,NULL,NULL,NULL,ts);CHKERRQ(ierr); 955 PetscFunctionReturn(0); 956 } 957 958 /*@C 959 TSGLEESetType - Set the type of GLEE scheme 960 961 Logically collective 962 963 Input Parameter: 964 + ts - timestepping context 965 - gleetype - type of GLEE-scheme 966 967 Level: intermediate 968 969 .seealso: TSGLEEGetType(), TSGLEE 970 @*/ 971 PetscErrorCode TSGLEESetType(TS ts,TSGLEEType gleetype) 972 { 973 PetscErrorCode ierr; 974 975 PetscFunctionBegin; 976 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 977 PetscValidCharPointer(gleetype,2); 978 ierr = PetscTryMethod(ts,"TSGLEESetType_C",(TS,TSGLEEType),(ts,gleetype));CHKERRQ(ierr); 979 PetscFunctionReturn(0); 980 } 981 982 /*@C 983 TSGLEEGetType - Get the type of GLEE scheme 984 985 Logically collective 986 987 Input Parameter: 988 . ts - timestepping context 989 990 Output Parameter: 991 . gleetype - type of GLEE-scheme 992 993 Level: intermediate 994 995 .seealso: TSGLEESetType() 996 @*/ 997 PetscErrorCode TSGLEEGetType(TS ts,TSGLEEType *gleetype) 998 { 999 PetscErrorCode ierr; 1000 1001 PetscFunctionBegin; 1002 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1003 ierr = PetscUseMethod(ts,"TSGLEEGetType_C",(TS,TSGLEEType*),(ts,gleetype));CHKERRQ(ierr); 1004 PetscFunctionReturn(0); 1005 } 1006 1007 PetscErrorCode TSGLEEGetType_GLEE(TS ts,TSGLEEType *gleetype) 1008 { 1009 TS_GLEE *glee = (TS_GLEE*)ts->data; 1010 PetscErrorCode ierr; 1011 1012 PetscFunctionBegin; 1013 if (!glee->tableau) { 1014 ierr = TSGLEESetType(ts,TSGLEEDefaultType);CHKERRQ(ierr); 1015 } 1016 *gleetype = glee->tableau->name; 1017 PetscFunctionReturn(0); 1018 } 1019 PetscErrorCode TSGLEESetType_GLEE(TS ts,TSGLEEType gleetype) 1020 { 1021 TS_GLEE *glee = (TS_GLEE*)ts->data; 1022 PetscErrorCode ierr; 1023 PetscBool match; 1024 GLEETableauLink link; 1025 1026 PetscFunctionBegin; 1027 if (glee->tableau) { 1028 ierr = PetscStrcmp(glee->tableau->name,gleetype,&match);CHKERRQ(ierr); 1029 if (match) PetscFunctionReturn(0); 1030 } 1031 for (link = GLEETableauList; link; link=link->next) { 1032 ierr = PetscStrcmp(link->tab.name,gleetype,&match);CHKERRQ(ierr); 1033 if (match) { 1034 ierr = TSReset_GLEE(ts);CHKERRQ(ierr); 1035 glee->tableau = &link->tab; 1036 PetscFunctionReturn(0); 1037 } 1038 } 1039 SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_UNKNOWN_TYPE,"Could not find '%s'",gleetype); 1040 PetscFunctionReturn(0); 1041 } 1042 1043 static PetscErrorCode TSGetStages_GLEE(TS ts,PetscInt *ns,Vec **Y) 1044 { 1045 TS_GLEE *glee = (TS_GLEE*)ts->data; 1046 1047 PetscFunctionBegin; 1048 if (ns) *ns = glee->tableau->s; 1049 if (Y) *Y = glee->YStage; 1050 PetscFunctionReturn(0); 1051 } 1052 1053 PetscErrorCode TSGetSolutionComponents_GLEE(TS ts,PetscInt *n,Vec *Y) 1054 { 1055 TS_GLEE *glee = (TS_GLEE*)ts->data; 1056 GLEETableau tab = glee->tableau; 1057 PetscErrorCode ierr; 1058 1059 PetscFunctionBegin; 1060 if (!Y) *n = tab->r; 1061 else { 1062 if ((*n >= 0) && (*n < tab->r)) { 1063 ierr = VecCopy(glee->Y[*n],*Y); CHKERRQ(ierr); 1064 } else SETERRQ3(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_OUTOFRANGE,"Second argument (%d) out of range[%d,%d].",*n,0,tab->r-1); 1065 } 1066 PetscFunctionReturn(0); 1067 } 1068 1069 PetscErrorCode TSGetAuxSolution_GLEE(TS ts,Vec *X) 1070 { 1071 TS_GLEE *glee = (TS_GLEE*)ts->data; 1072 GLEETableau tab = glee->tableau; 1073 PetscReal *F = tab->Fembed; 1074 PetscInt r = tab->r; 1075 Vec *Y = glee->Y; 1076 PetscScalar *wr = glee->rwork; 1077 PetscInt i; 1078 PetscErrorCode ierr; 1079 1080 PetscFunctionBegin; 1081 ierr = VecZeroEntries(*X);CHKERRQ(ierr); 1082 for (i=0; i<r; i++) wr[i] = F[i]; 1083 ierr = VecMAXPY((*X),r,wr,Y);CHKERRQ(ierr); 1084 PetscFunctionReturn(0); 1085 } 1086 1087 PetscErrorCode TSGetTimeError_GLEE(TS ts,PetscInt n,Vec *X) 1088 { 1089 TS_GLEE *glee = (TS_GLEE*)ts->data; 1090 GLEETableau tab = glee->tableau; 1091 PetscReal *F = tab->Ferror; 1092 PetscInt r = tab->r; 1093 Vec *Y = glee->Y; 1094 PetscScalar *wr = glee->rwork; 1095 PetscInt i; 1096 PetscErrorCode ierr; 1097 1098 PetscFunctionBegin; 1099 ierr = VecZeroEntries(*X);CHKERRQ(ierr); 1100 if(n==0){ 1101 for (i=0; i<r; i++) wr[i] = F[i]; 1102 ierr = VecMAXPY((*X),r,wr,Y);CHKERRQ(ierr); 1103 } else if(n==-1) { 1104 *X=glee->yGErr; 1105 } 1106 PetscFunctionReturn(0); 1107 } 1108 1109 PetscErrorCode TSSetTimeError_GLEE(TS ts,Vec X) 1110 { 1111 TS_GLEE *glee = (TS_GLEE*)ts->data; 1112 GLEETableau tab = glee->tableau; 1113 PetscReal *S = tab->Serror; 1114 PetscInt r = tab->r,i; 1115 Vec *Y = glee->Y; 1116 PetscErrorCode ierr; 1117 1118 PetscFunctionBegin; 1119 if (r != 2) SETERRQ2(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"TSSetTimeError_GLEE not supported for '%s' with r=%D.",tab->name,tab->r); 1120 for (i=1; i<r; i++) { 1121 ierr = VecCopy(ts->vec_sol,Y[i]); CHKERRQ(ierr); 1122 ierr = VecAXPBY(Y[i],S[0],S[1],X); CHKERRQ(ierr); 1123 ierr = VecCopy(X,glee->yGErr); CHKERRQ(ierr); 1124 } 1125 PetscFunctionReturn(0); 1126 } 1127 1128 static PetscErrorCode TSDestroy_GLEE(TS ts) 1129 { 1130 PetscErrorCode ierr; 1131 1132 PetscFunctionBegin; 1133 ierr = TSReset_GLEE(ts);CHKERRQ(ierr); 1134 if (ts->dm) { 1135 ierr = DMCoarsenHookRemove(ts->dm,DMCoarsenHook_TSGLEE,DMRestrictHook_TSGLEE,ts);CHKERRQ(ierr); 1136 ierr = DMSubDomainHookRemove(ts->dm,DMSubDomainHook_TSGLEE,DMSubDomainRestrictHook_TSGLEE,ts);CHKERRQ(ierr); 1137 } 1138 ierr = PetscFree(ts->data);CHKERRQ(ierr); 1139 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSGLEEGetType_C",NULL);CHKERRQ(ierr); 1140 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSGLEESetType_C",NULL);CHKERRQ(ierr); 1141 PetscFunctionReturn(0); 1142 } 1143 1144 /* ------------------------------------------------------------ */ 1145 /*MC 1146 TSGLEE - ODE and DAE solver using General Linear with Error Estimation schemes 1147 1148 The user should provide the right hand side of the equation 1149 using TSSetRHSFunction(). 1150 1151 Notes: 1152 The default is TSGLEE35, it can be changed with TSGLEESetType() or -ts_glee_type 1153 1154 Level: beginner 1155 1156 .seealso: TSCreate(), TS, TSSetType(), TSGLEESetType(), TSGLEEGetType(), 1157 TSGLEE23, TTSGLEE24, TSGLEE35, TSGLEE25I, TSGLEEEXRK2A, 1158 TSGLEERK32G1, TSGLEERK285EX, TSGLEEType, TSGLEERegister() 1159 1160 M*/ 1161 PETSC_EXTERN PetscErrorCode TSCreate_GLEE(TS ts) 1162 { 1163 TS_GLEE *th; 1164 PetscErrorCode ierr; 1165 1166 PetscFunctionBegin; 1167 #if !defined(PETSC_USE_DYNAMIC_LIBRARIES) 1168 ierr = TSGLEEInitializePackage();CHKERRQ(ierr); 1169 #endif 1170 1171 ts->ops->reset = TSReset_GLEE; 1172 ts->ops->destroy = TSDestroy_GLEE; 1173 ts->ops->view = TSView_GLEE; 1174 ts->ops->load = TSLoad_GLEE; 1175 ts->ops->setup = TSSetUp_GLEE; 1176 ts->ops->step = TSStep_GLEE; 1177 ts->ops->interpolate = TSInterpolate_GLEE; 1178 ts->ops->evaluatestep = TSEvaluateStep_GLEE; 1179 ts->ops->setfromoptions = TSSetFromOptions_GLEE; 1180 ts->ops->getstages = TSGetStages_GLEE; 1181 ts->ops->snesfunction = SNESTSFormFunction_GLEE; 1182 ts->ops->snesjacobian = SNESTSFormJacobian_GLEE; 1183 ts->ops->getsolutioncomponents = TSGetSolutionComponents_GLEE; 1184 ts->ops->getauxsolution = TSGetAuxSolution_GLEE; 1185 ts->ops->gettimeerror = TSGetTimeError_GLEE; 1186 ts->ops->settimeerror = TSSetTimeError_GLEE; 1187 ts->ops->startingmethod = TSStartingMethod_GLEE; 1188 ts->default_adapt_type = TSADAPTGLEE; 1189 1190 ts->usessnes = PETSC_TRUE; 1191 1192 ierr = PetscNewLog(ts,&th);CHKERRQ(ierr); 1193 ts->data = (void*)th; 1194 1195 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSGLEEGetType_C",TSGLEEGetType_GLEE);CHKERRQ(ierr); 1196 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSGLEESetType_C",TSGLEESetType_GLEE);CHKERRQ(ierr); 1197 PetscFunctionReturn(0); 1198 } 1199