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