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