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