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