1 /* 2 Code for time stepping with the General Linear with Error Estimation method 3 4 Notes: 5 The general system is written as 6 7 Udot = F(t,U) 8 9 */ 10 #include <petsc/private/tsimpl.h> /*I "petscts.h" I*/ 11 #include <petscdm.h> 12 13 static PetscBool cited = PETSC_FALSE; 14 static const char citation[] = 15 "@ARTICLE{Constantinescu_TR2016b,\n" 16 " author = {Constantinescu, E.M.},\n" 17 " title = {Estimating Global Errors in Time Stepping},\n" 18 " journal = {ArXiv e-prints},\n" 19 " year = 2016,\n" 20 " adsurl = {http://adsabs.harvard.edu/abs/2015arXiv150305166C}\n}\n"; 21 22 static TSGLEEType TSGLEEDefaultType = TSGLEE35; 23 static PetscBool TSGLEERegisterAllCalled; 24 static PetscBool TSGLEEPackageInitialized; 25 static PetscInt explicit_stage_time_id; 26 27 typedef struct _GLEETableau *GLEETableau; 28 struct _GLEETableau { 29 char *name; 30 PetscInt order; /* Classical approximation order of the method i*/ 31 PetscInt s; /* Number of stages */ 32 PetscInt r; /* Number of steps */ 33 PetscReal gamma; /* LTE ratio */ 34 PetscReal *A,*B,*U,*V,*S,*F,*c; /* Tableau */ 35 PetscReal *Fembed; /* Embedded final method coefficients */ 36 PetscReal *Ferror; /* Coefficients for computing error */ 37 PetscReal *Serror; /* Coefficients for initializing the error */ 38 PetscInt pinterp; /* Interpolation order */ 39 PetscReal *binterp; /* Interpolation coefficients */ 40 PetscReal ccfl; /* Placeholder for CFL coefficient relative to forward Euler */ 41 }; 42 typedef struct _GLEETableauLink *GLEETableauLink; 43 struct _GLEETableauLink { 44 struct _GLEETableau tab; 45 GLEETableauLink next; 46 }; 47 static GLEETableauLink GLEETableauList; 48 49 typedef struct { 50 GLEETableau tableau; 51 Vec *Y; /* Solution vector (along with auxiliary solution y~ or eps) */ 52 Vec *X; /* Temporary solution vector */ 53 Vec *YStage; /* Stage values */ 54 Vec *YdotStage; /* Stage right hand side */ 55 Vec W; /* Right-hand-side for implicit stage solve */ 56 Vec Ydot; /* Work vector holding Ydot during residual evaluation */ 57 Vec yGErr; /* Vector holding the global error after a step is completed */ 58 PetscScalar *swork; /* Scalar work (size of the number of stages)*/ 59 PetscScalar *rwork; /* Scalar work (size of the number of steps)*/ 60 PetscReal scoeff; /* shift = scoeff/dt */ 61 PetscReal stage_time; 62 TSStepStatus status; 63 } TS_GLEE; 64 65 /*MC 66 TSGLEE23 - Second order three stage GLEE method 67 68 This method has three stages. 69 s = 3, r = 2 70 71 Level: advanced 72 73 .seealso: TSGLEE 74 M*/ 75 /*MC 76 TSGLEE24 - Second order four stage GLEE method 77 78 This method has four stages. 79 s = 4, r = 2 80 81 Level: advanced 82 83 .seealso: TSGLEE 84 M*/ 85 /*MC 86 TSGLEE25i - Second order five stage GLEE method 87 88 This method has five stages. 89 s = 5, r = 2 90 91 Level: advanced 92 93 .seealso: TSGLEE 94 M*/ 95 /*MC 96 TSGLEE35 - Third order five stage GLEE method 97 98 This method has five stages. 99 s = 5, r = 2 100 101 Level: advanced 102 103 .seealso: TSGLEE 104 M*/ 105 /*MC 106 TSGLEEEXRK2A - Second order six stage GLEE method 107 108 This method has six stages. 109 s = 6, r = 2 110 111 Level: advanced 112 113 .seealso: TSGLEE 114 M*/ 115 /*MC 116 TSGLEERK32G1 - Third order eight stage GLEE method 117 118 This method has eight stages. 119 s = 8, r = 2 120 121 Level: advanced 122 123 .seealso: TSGLEE 124 M*/ 125 /*MC 126 TSGLEERK285EX - Second order nine stage GLEE method 127 128 This method has nine stages. 129 s = 9, r = 2 130 131 Level: advanced 132 133 .seealso: TSGLEE 134 M*/ 135 136 /*@C 137 TSGLEERegisterAll - Registers all of the General Linear with Error Estimation methods in TSGLEE 138 139 Not Collective, but should be called by all processes which will need the schemes to be registered 140 141 Level: advanced 142 143 .seealso: TSGLEERegisterDestroy() 144 @*/ 145 PetscErrorCode TSGLEERegisterAll(void) 146 { 147 PetscErrorCode ierr; 148 149 PetscFunctionBegin; 150 if (TSGLEERegisterAllCalled) PetscFunctionReturn(0); 151 TSGLEERegisterAllCalled = PETSC_TRUE; 152 153 { 154 #define GAMMA 0.5 155 /* y-eps form */ 156 const PetscInt 157 p = 1, 158 s = 3, 159 r = 2; 160 const PetscReal 161 A[3][3] = {{1.0,0,0},{0,0.5,0},{0,0.5,0.5}}, 162 B[2][3] = {{1.0,0,0},{-2.0,1.0,1.0}}, 163 U[3][2] = {{1.0,0},{1.0,0.5},{1.0,0.5}}, 164 V[2][2] = {{1,0},{0,1}}, 165 S[2] = {1,0}, 166 F[2] = {1,0}, 167 Fembed[2] = {1,1-GAMMA}, 168 Ferror[2] = {0,1}, 169 Serror[2] = {1,0}; 170 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); 171 } 172 { 173 #undef GAMMA 174 #define GAMMA 0.0 175 /* y-eps form */ 176 const PetscInt 177 p = 2, 178 s = 3, 179 r = 2; 180 const PetscReal 181 A[3][3] = {{0,0,0},{1,0,0},{0.25,0.25,0}}, 182 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}}, 183 U[3][2] = {{1,0},{1,10},{1,-1}}, 184 V[2][2] = {{1,0},{0,1}}, 185 S[2] = {1,0}, 186 F[2] = {1,0}, 187 Fembed[2] = {1,1-GAMMA}, 188 Ferror[2] = {0,1}, 189 Serror[2] = {1,0}; 190 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); 191 } 192 { 193 #undef GAMMA 194 #define GAMMA 0.0 195 /* y-y~ form */ 196 const PetscInt 197 p = 2, 198 s = 4, 199 r = 2; 200 const PetscReal 201 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}}, 202 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}}, 203 U[4][2] = {{0,1},{75.0/58.0,-17.0/58.0},{0,1},{0,1}}, 204 V[2][2] = {{1,0},{0,1}}, 205 S[2] = {1,1}, 206 F[2] = {1,0}, 207 Fembed[2] = {0,1}, 208 Ferror[2] = {-1.0/(1.0-GAMMA),1.0/(1.0-GAMMA)}, 209 Serror[2] = {1.0-GAMMA,1.0}; 210 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); 211 } 212 { 213 #undef GAMMA 214 #define GAMMA 0.0 215 /* y-y~ form */ 216 const PetscInt 217 p = 2, 218 s = 5, 219 r = 2; 220 const PetscReal 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 #undef GAMMA 243 #define GAMMA 0.0 244 /* y-y~ form */ 245 const PetscInt 246 p = 3, 247 s = 5, 248 r = 2; 249 const PetscReal 250 A[5][5] = {{0,0,0,0,0}, 251 {- 2169604947363702313.0 / 24313474998937147335.0,0,0,0,0}, 252 {46526746497697123895.0 / 94116917485856474137.0,-10297879244026594958.0 / 49199457603717988219.0,0,0,0}, 253 {23364788935845982499.0 / 87425311444725389446.0,-79205144337496116638.0 / 148994349441340815519.0,40051189859317443782.0 / 36487615018004984309.0,0,0}, 254 {42089522664062539205.0 / 124911313006412840286.0,-15074384760342762939.0 / 137927286865289746282.0,-62274678522253371016.0 / 125918573676298591413.0,13755475729852471739.0 / 79257927066651693390.0,0}}, 255 B[2][5] = {{61546696837458703723.0 / 56982519523786160813.0,-55810892792806293355.0 / 206957624151308356511.0,24061048952676379087.0 / 158739347956038723465.0,3577972206874351339.0 / 7599733370677197135.0,-59449832954780563947.0 / 137360038685338563670.0}, 256 {- 9738262186984159168.0 / 99299082461487742983.0,-32797097931948613195.0 / 61521565616362163366.0,42895514606418420631.0 / 71714201188501437336.0,22608567633166065068.0 / 55371917805607957003.0,94655809487476459565.0 / 151517167160302729021.0}}, 257 U[5][2] = {{70820309139834661559.0 / 80863923579509469826.0,10043614439674808267.0 / 80863923579509469826.0}, 258 {161694774978034105510.0 / 106187653640211060371.0,-55507121337823045139.0 / 106187653640211060371.0}, 259 {78486094644566264568.0 / 88171030896733822981.0,9684936252167558413.0 / 88171030896733822981.0}, 260 {65394922146334854435.0 / 84570853840405479554.0,19175931694070625119.0 / 84570853840405479554.0}, 261 {8607282770183754108.0 / 108658046436496925911.0,100050763666313171803.0 / 108658046436496925911.0}}, 262 V[2][2] = {{1,0},{0,1}}, 263 S[2] = {1,1}, 264 F[2] = {1,0}, 265 Fembed[2] = {0,1}, 266 Ferror[2] = {-1.0/(1.0-GAMMA),1.0/(1.0-GAMMA)}, 267 Serror[2] = {1.0-GAMMA,1.0}; 268 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); 269 } 270 { 271 #undef GAMMA 272 #define GAMMA 0.25 273 /* y-eps form */ 274 const PetscInt 275 p = 2, 276 s = 6, 277 r = 2; 278 const PetscReal 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 #undef GAMMA 298 #define GAMMA 0.0 299 /* y-eps form */ 300 const PetscInt 301 p = 3, 302 s = 8, 303 r = 2; 304 const PetscReal 305 A[8][8] = {{0,0,0,0,0,0,0,0}, 306 {0.5,0,0,0,0,0,0,0}, 307 {-1,2,0,0,0,0,0,0}, 308 {1.0/6.0,2.0/3.0,1.0/6.0,0,0,0,0,0}, 309 {0,0,0,0,0,0,0,0}, 310 {-7.0/24.0,1.0/3.0,1.0/12.0,-1.0/8.0,0.5,0,0,0}, 311 {7.0/6.0,-4.0/3.0,-1.0/3.0,0.5,-1.0,2.0,0,0}, 312 {0,0,0,0,1.0/6.0,2.0/3.0,1.0/6.0,0}}, 313 B[2][8] = {{1.0/6.0,2.0/3.0,1.0/6.0,0,0,0,0,0}, 314 {-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}}, 315 U[8][2] = {{1,0},{1,0},{1,0},{1,0},{1,1},{1,1},{1,1},{1,1}}, 316 V[2][2] = {{1,0},{0,1}}, 317 S[2] = {1,0}, 318 F[2] = {1,0}, 319 Fembed[2] = {1,1-GAMMA}, 320 Ferror[2] = {0,1}, 321 Serror[2] = {1,0}; 322 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); 323 } 324 { 325 #undef GAMMA 326 #define GAMMA 0.25 327 /* y-eps form */ 328 const PetscInt 329 p = 2, 330 s = 9, 331 r = 2; 332 const PetscReal 333 A[9][9] = {{0,0,0,0,0,0,0,0,0}, 334 {0.585786437626904966,0,0,0,0,0,0,0,0}, 335 {0.149999999999999994,0.849999999999999978,0,0,0,0,0,0,0}, 336 {0,0,0,0,0,0,0,0,0}, 337 {0,0,0,0.292893218813452483,0,0,0,0,0}, 338 {0,0,0,0.0749999999999999972,0.424999999999999989,0,0,0,0}, 339 {0,0,0,0.176776695296636893,0.176776695296636893,0.146446609406726241,0,0,0}, 340 {0,0,0,0.176776695296636893,0.176776695296636893,0.146446609406726241,0.292893218813452483,0,0}, 341 {0,0,0,0.176776695296636893,0.176776695296636893,0.146446609406726241,0.0749999999999999972,0.424999999999999989,0}}, 342 B[2][9] = {{0.353553390593273786,0.353553390593273786,0.292893218813452483,0,0,0,0,0,0}, 343 {-0.471404520791031678,-0.471404520791031678,-0.390524291751269959,0.235702260395515839,0.235702260395515839,0.195262145875634979,0.235702260395515839,0.235702260395515839,0.195262145875634979}}, 344 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}}, 345 V[2][2] = {{1,0},{0,1}}, 346 S[2] = {1,0}, 347 F[2] = {1,0}, 348 Fembed[2] = {1,1-GAMMA}, 349 Ferror[2] = {0,1}, 350 Serror[2] = {1,0}; 351 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); 352 } 353 354 PetscFunctionReturn(0); 355 } 356 357 /*@C 358 TSGLEERegisterDestroy - Frees the list of schemes that were registered by TSGLEERegister(). 359 360 Not Collective 361 362 Level: advanced 363 364 .seealso: TSGLEERegister(), TSGLEERegisterAll() 365 @*/ 366 PetscErrorCode TSGLEERegisterDestroy(void) 367 { 368 PetscErrorCode ierr; 369 GLEETableauLink link; 370 371 PetscFunctionBegin; 372 while ((link = GLEETableauList)) { 373 GLEETableau t = &link->tab; 374 GLEETableauList = link->next; 375 ierr = PetscFree5(t->A,t->B,t->U,t->V,t->c);CHKERRQ(ierr); 376 ierr = PetscFree2(t->S,t->F); CHKERRQ(ierr); 377 ierr = PetscFree (t->Fembed); CHKERRQ(ierr); 378 ierr = PetscFree (t->Ferror); CHKERRQ(ierr); 379 ierr = PetscFree (t->Serror); CHKERRQ(ierr); 380 ierr = PetscFree (t->binterp); CHKERRQ(ierr); 381 ierr = PetscFree (t->name); CHKERRQ(ierr); 382 ierr = PetscFree (link); CHKERRQ(ierr); 383 } 384 TSGLEERegisterAllCalled = PETSC_FALSE; 385 PetscFunctionReturn(0); 386 } 387 388 /*@C 389 TSGLEEInitializePackage - This function initializes everything in the TSGLEE package. It is called 390 from TSInitializePackage(). 391 392 Level: developer 393 394 .seealso: PetscInitialize() 395 @*/ 396 PetscErrorCode TSGLEEInitializePackage(void) 397 { 398 PetscErrorCode ierr; 399 400 PetscFunctionBegin; 401 if (TSGLEEPackageInitialized) PetscFunctionReturn(0); 402 TSGLEEPackageInitialized = PETSC_TRUE; 403 ierr = TSGLEERegisterAll();CHKERRQ(ierr); 404 ierr = PetscObjectComposedDataRegister(&explicit_stage_time_id);CHKERRQ(ierr); 405 ierr = PetscRegisterFinalize(TSGLEEFinalizePackage);CHKERRQ(ierr); 406 PetscFunctionReturn(0); 407 } 408 409 /*@C 410 TSGLEEFinalizePackage - This function destroys everything in the TSGLEE package. It is 411 called from PetscFinalize(). 412 413 Level: developer 414 415 .seealso: PetscFinalize() 416 @*/ 417 PetscErrorCode TSGLEEFinalizePackage(void) 418 { 419 PetscErrorCode ierr; 420 421 PetscFunctionBegin; 422 TSGLEEPackageInitialized = PETSC_FALSE; 423 ierr = TSGLEERegisterDestroy();CHKERRQ(ierr); 424 PetscFunctionReturn(0); 425 } 426 427 /*@C 428 TSGLEERegister - register an GLEE scheme by providing the entries in the Butcher tableau 429 430 Not Collective, but the same schemes should be registered on all processes on which they will be used 431 432 Input Parameters: 433 + name - identifier for method 434 . order - order of method 435 . s - number of stages 436 . r - number of steps 437 . gamma - LTE ratio 438 . A - stage coefficients (dimension s*s, row-major) 439 . B - step completion coefficients (dimension r*s, row-major) 440 . U - method coefficients (dimension s*r, row-major) 441 . V - method coefficients (dimension r*r, row-major) 442 . S - starting coefficients 443 . F - finishing coefficients 444 . c - abscissa (dimension s; NULL to use row sums of A) 445 . Fembed - step completion coefficients for embedded method 446 . Ferror - error computation coefficients 447 . Serror - error initialization coefficients 448 . pinterp - order of interpolation (0 if unavailable) 449 - binterp - array of interpolation coefficients (NULL if unavailable) 450 451 Notes: 452 Several GLEE methods are provided, this function is only needed to create new methods. 453 454 Level: advanced 455 456 .seealso: TSGLEE 457 @*/ 458 PetscErrorCode TSGLEERegister(TSGLEEType name,PetscInt order,PetscInt s, PetscInt r, 459 PetscReal gamma, 460 const PetscReal A[],const PetscReal B[], 461 const PetscReal U[],const PetscReal V[], 462 const PetscReal S[],const PetscReal F[], 463 const PetscReal c[], 464 const PetscReal Fembed[],const PetscReal Ferror[], 465 const PetscReal Serror[], 466 PetscInt pinterp, const PetscReal binterp[]) 467 { 468 PetscErrorCode ierr; 469 GLEETableauLink link; 470 GLEETableau t; 471 PetscInt i,j; 472 473 PetscFunctionBegin; 474 ierr = TSGLEEInitializePackage();CHKERRQ(ierr); 475 ierr = PetscNew(&link);CHKERRQ(ierr); 476 t = &link->tab; 477 ierr = PetscStrallocpy(name,&t->name);CHKERRQ(ierr); 478 t->order = order; 479 t->s = s; 480 t->r = r; 481 t->gamma = gamma; 482 ierr = PetscMalloc5(s*s,&t->A,r*r,&t->V,s,&t->c,r*s,&t->B,s*r,&t->U);CHKERRQ(ierr); 483 ierr = PetscMalloc2(r,&t->S,r,&t->F);CHKERRQ(ierr); 484 ierr = PetscArraycpy(t->A,A,s*s);CHKERRQ(ierr); 485 ierr = PetscArraycpy(t->B,B,r*s);CHKERRQ(ierr); 486 ierr = PetscArraycpy(t->U,U,s*r);CHKERRQ(ierr); 487 ierr = PetscArraycpy(t->V,V,r*r);CHKERRQ(ierr); 488 ierr = PetscArraycpy(t->S,S,r);CHKERRQ(ierr); 489 ierr = PetscArraycpy(t->F,F,r);CHKERRQ(ierr); 490 if (c) { 491 ierr = PetscArraycpy(t->c,c,s);CHKERRQ(ierr); 492 } else { 493 for (i=0; i<s; i++) for (j=0,t->c[i]=0; j<s; j++) t->c[i] += A[i*s+j]; 494 } 495 ierr = PetscMalloc1(r,&t->Fembed);CHKERRQ(ierr); 496 ierr = PetscMalloc1(r,&t->Ferror);CHKERRQ(ierr); 497 ierr = PetscMalloc1(r,&t->Serror);CHKERRQ(ierr); 498 ierr = PetscArraycpy(t->Fembed,Fembed,r);CHKERRQ(ierr); 499 ierr = PetscArraycpy(t->Ferror,Ferror,r);CHKERRQ(ierr); 500 ierr = PetscArraycpy(t->Serror,Serror,r);CHKERRQ(ierr); 501 t->pinterp = pinterp; 502 ierr = PetscMalloc1(s*pinterp,&t->binterp);CHKERRQ(ierr); 503 ierr = PetscArraycpy(t->binterp,binterp,s*pinterp);CHKERRQ(ierr); 504 505 link->next = GLEETableauList; 506 GLEETableauList = link; 507 PetscFunctionReturn(0); 508 } 509 510 static PetscErrorCode TSEvaluateStep_GLEE(TS ts,PetscInt order,Vec X,PetscBool *done) 511 { 512 TS_GLEE *glee = (TS_GLEE*) ts->data; 513 GLEETableau tab = glee->tableau; 514 PetscReal h, *B = tab->B, *V = tab->V, 515 *F = tab->F, 516 *Fembed = tab->Fembed; 517 PetscInt s = tab->s, r = tab->r, i, j; 518 Vec *Y = glee->Y, *YdotStage = glee->YdotStage; 519 PetscScalar *ws = glee->swork, *wr = glee->rwork; 520 PetscErrorCode ierr; 521 522 PetscFunctionBegin; 523 524 switch (glee->status) { 525 case TS_STEP_INCOMPLETE: 526 case TS_STEP_PENDING: 527 h = ts->time_step; break; 528 case TS_STEP_COMPLETE: 529 h = ts->ptime - ts->ptime_prev; break; 530 default: SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_PLIB,"Invalid TSStepStatus"); 531 } 532 533 if (order == tab->order) { 534 535 /* Note: Irrespective of whether status is TS_STEP_INCOMPLETE 536 or TS_STEP_COMPLETE, glee->X has the solution at the 537 beginning of the time step. So no need to roll-back. 538 */ 539 if (glee->status == TS_STEP_INCOMPLETE) { 540 for (i=0; i<r; i++) { 541 ierr = VecZeroEntries(Y[i]);CHKERRQ(ierr); 542 for (j=0; j<r; j++) wr[j] = V[i*r+j]; 543 ierr = VecMAXPY(Y[i],r,wr,glee->X);CHKERRQ(ierr); 544 for (j=0; j<s; j++) ws[j] = h*B[i*s+j]; 545 ierr = VecMAXPY(Y[i],s,ws,YdotStage);CHKERRQ(ierr); 546 } 547 ierr = VecZeroEntries(X);CHKERRQ(ierr); 548 for (j=0; j<r; j++) wr[j] = F[j]; 549 ierr = VecMAXPY(X,r,wr,Y);CHKERRQ(ierr); 550 } else {ierr = VecCopy(ts->vec_sol,X);CHKERRQ(ierr);} 551 PetscFunctionReturn(0); 552 553 } else if (order == tab->order-1) { 554 555 /* Complete with the embedded method (Fembed) */ 556 for (i=0; i<r; i++) { 557 ierr = VecZeroEntries(Y[i]);CHKERRQ(ierr); 558 for (j=0; j<r; j++) wr[j] = V[i*r+j]; 559 ierr = VecMAXPY(Y[i],r,wr,glee->X);CHKERRQ(ierr); 560 for (j=0; j<s; j++) ws[j] = h*B[i*s+j]; 561 ierr = VecMAXPY(Y[i],s,ws,YdotStage);CHKERRQ(ierr); 562 } 563 ierr = VecZeroEntries(X);CHKERRQ(ierr); 564 for (j=0; j<r; j++) wr[j] = Fembed[j]; 565 ierr = VecMAXPY(X,r,wr,Y);CHKERRQ(ierr); 566 567 if (done) *done = PETSC_TRUE; 568 PetscFunctionReturn(0); 569 } 570 if (done) *done = PETSC_FALSE; 571 else SETERRQ3(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"GLEE '%s' of order %D cannot evaluate step at order %D",tab->name,tab->order,order); 572 PetscFunctionReturn(0); 573 } 574 575 static PetscErrorCode TSStep_GLEE(TS ts) 576 { 577 TS_GLEE *glee = (TS_GLEE*)ts->data; 578 GLEETableau tab = glee->tableau; 579 const PetscInt s = tab->s, r = tab->r; 580 PetscReal *A = tab->A, *U = tab->U, 581 *F = tab->F, 582 *c = tab->c; 583 Vec *Y = glee->Y, *X = glee->X, 584 *YStage = glee->YStage, 585 *YdotStage = glee->YdotStage, 586 W = glee->W; 587 SNES snes; 588 PetscScalar *ws = glee->swork, *wr = glee->rwork; 589 TSAdapt adapt; 590 PetscInt i,j,reject,next_scheme,its,lits; 591 PetscReal next_time_step; 592 PetscReal t; 593 PetscBool accept; 594 PetscErrorCode ierr; 595 596 PetscFunctionBegin; 597 ierr = PetscCitationsRegister(citation,&cited);CHKERRQ(ierr); 598 599 for (i=0; i<r; i++) { ierr = VecCopy(Y[i],X[i]);CHKERRQ(ierr); } 600 601 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 602 next_time_step = ts->time_step; 603 t = ts->ptime; 604 accept = PETSC_TRUE; 605 glee->status = TS_STEP_INCOMPLETE; 606 607 for (reject=0; reject<ts->max_reject && !ts->reason; reject++,ts->reject++) { 608 609 PetscReal h = ts->time_step; 610 ierr = TSPreStep(ts);CHKERRQ(ierr); 611 612 for (i=0; i<s; i++) { 613 614 glee->stage_time = t + h*c[i]; 615 ierr = TSPreStage(ts,glee->stage_time);CHKERRQ(ierr); 616 617 if (A[i*s+i] == 0) { /* Explicit stage */ 618 ierr = VecZeroEntries(YStage[i]);CHKERRQ(ierr); 619 for (j=0; j<r; j++) wr[j] = U[i*r+j]; 620 ierr = VecMAXPY(YStage[i],r,wr,X);CHKERRQ(ierr); 621 for (j=0; j<i; j++) ws[j] = h*A[i*s+j]; 622 ierr = VecMAXPY(YStage[i],i,ws,YdotStage);CHKERRQ(ierr); 623 } else { /* Implicit stage */ 624 glee->scoeff = 1.0/A[i*s+i]; 625 /* compute right-hand-side */ 626 ierr = VecZeroEntries(W);CHKERRQ(ierr); 627 for (j=0; j<r; j++) wr[j] = U[i*r+j]; 628 ierr = VecMAXPY(W,r,wr,X);CHKERRQ(ierr); 629 for (j=0; j<i; j++) ws[j] = h*A[i*s+j]; 630 ierr = VecMAXPY(W,i,ws,YdotStage);CHKERRQ(ierr); 631 ierr = VecScale(W,glee->scoeff/h);CHKERRQ(ierr); 632 /* set initial guess */ 633 ierr = VecCopy(i>0 ? YStage[i-1] : ts->vec_sol,YStage[i]);CHKERRQ(ierr); 634 /* solve for this stage */ 635 ierr = SNESSolve(snes,W,YStage[i]);CHKERRQ(ierr); 636 ierr = SNESGetIterationNumber(snes,&its);CHKERRQ(ierr); 637 ierr = SNESGetLinearSolveIterations(snes,&lits);CHKERRQ(ierr); 638 ts->snes_its += its; ts->ksp_its += lits; 639 } 640 ierr = TSGetAdapt(ts,&adapt);CHKERRQ(ierr); 641 ierr = TSAdaptCheckStage(adapt,ts,glee->stage_time,YStage[i],&accept);CHKERRQ(ierr); 642 if (!accept) goto reject_step; 643 ierr = TSPostStage(ts,glee->stage_time,i,YStage);CHKERRQ(ierr); 644 ierr = TSComputeRHSFunction(ts,t+h*c[i],YStage[i],YdotStage[i]);CHKERRQ(ierr); 645 } 646 ierr = TSEvaluateStep(ts,tab->order,ts->vec_sol,NULL);CHKERRQ(ierr); 647 glee->status = TS_STEP_PENDING; 648 649 /* Register only the current method as a candidate because we're not supporting multiple candidates yet. */ 650 ierr = TSGetAdapt(ts,&adapt);CHKERRQ(ierr); 651 ierr = TSAdaptCandidatesClear(adapt);CHKERRQ(ierr); 652 ierr = TSAdaptCandidateAdd(adapt,tab->name,tab->order,1,tab->ccfl,(PetscReal)tab->s,PETSC_TRUE);CHKERRQ(ierr); 653 ierr = TSAdaptChoose(adapt,ts,ts->time_step,&next_scheme,&next_time_step,&accept);CHKERRQ(ierr); 654 if (accept) { 655 /* ignore next_scheme for now */ 656 ts->ptime += ts->time_step; 657 ts->time_step = next_time_step; 658 glee->status = TS_STEP_COMPLETE; 659 /* compute and store the global error */ 660 /* Note: this is not needed if TSAdaptGLEE is not used */ 661 ierr = TSGetTimeError(ts,0,&(glee->yGErr));CHKERRQ(ierr); 662 ierr = PetscObjectComposedDataSetReal((PetscObject)ts->vec_sol,explicit_stage_time_id,ts->ptime);CHKERRQ(ierr); 663 break; 664 } else { /* Roll back the current step */ 665 for (j=0; j<r; j++) wr[j] = F[j]; 666 ierr = VecMAXPY(ts->vec_sol,r,wr,X);CHKERRQ(ierr); 667 ts->time_step = next_time_step; 668 glee->status = TS_STEP_INCOMPLETE; 669 } 670 reject_step: continue; 671 } 672 if (glee->status != TS_STEP_COMPLETE && !ts->reason) ts->reason = TS_DIVERGED_STEP_REJECTED; 673 PetscFunctionReturn(0); 674 } 675 676 static PetscErrorCode TSInterpolate_GLEE(TS ts,PetscReal itime,Vec X) 677 { 678 TS_GLEE *glee = (TS_GLEE*)ts->data; 679 PetscInt s=glee->tableau->s, pinterp=glee->tableau->pinterp,i,j; 680 PetscReal h,tt,t; 681 PetscScalar *b; 682 const PetscReal *B = glee->tableau->binterp; 683 PetscErrorCode ierr; 684 685 PetscFunctionBegin; 686 if (!B) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"TSGLEE %s does not have an interpolation formula",glee->tableau->name); 687 switch (glee->status) { 688 case TS_STEP_INCOMPLETE: 689 case TS_STEP_PENDING: 690 h = ts->time_step; 691 t = (itime - ts->ptime)/h; 692 break; 693 case TS_STEP_COMPLETE: 694 h = ts->ptime - ts->ptime_prev; 695 t = (itime - ts->ptime)/h + 1; /* In the interval [0,1] */ 696 break; 697 default: SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_PLIB,"Invalid TSStepStatus"); 698 } 699 ierr = PetscMalloc1(s,&b);CHKERRQ(ierr); 700 for (i=0; i<s; i++) b[i] = 0; 701 for (j=0,tt=t; j<pinterp; j++,tt*=t) { 702 for (i=0; i<s; i++) { 703 b[i] += h * B[i*pinterp+j] * tt; 704 } 705 } 706 ierr = VecCopy(glee->YStage[0],X);CHKERRQ(ierr); 707 ierr = VecMAXPY(X,s,b,glee->YdotStage);CHKERRQ(ierr); 708 ierr = PetscFree(b);CHKERRQ(ierr); 709 PetscFunctionReturn(0); 710 } 711 712 /*------------------------------------------------------------*/ 713 static PetscErrorCode TSReset_GLEE(TS ts) 714 { 715 TS_GLEE *glee = (TS_GLEE*)ts->data; 716 PetscInt s, r; 717 PetscErrorCode ierr; 718 719 PetscFunctionBegin; 720 if (!glee->tableau) PetscFunctionReturn(0); 721 s = glee->tableau->s; 722 r = glee->tableau->r; 723 ierr = VecDestroyVecs(r,&glee->Y);CHKERRQ(ierr); 724 ierr = VecDestroyVecs(r,&glee->X);CHKERRQ(ierr); 725 ierr = VecDestroyVecs(s,&glee->YStage);CHKERRQ(ierr); 726 ierr = VecDestroyVecs(s,&glee->YdotStage);CHKERRQ(ierr); 727 ierr = VecDestroy(&glee->Ydot);CHKERRQ(ierr); 728 ierr = VecDestroy(&glee->yGErr);CHKERRQ(ierr); 729 ierr = VecDestroy(&glee->W);CHKERRQ(ierr); 730 ierr = PetscFree2(glee->swork,glee->rwork);CHKERRQ(ierr); 731 PetscFunctionReturn(0); 732 } 733 734 static PetscErrorCode TSGLEEGetVecs(TS ts,DM dm,Vec *Ydot) 735 { 736 TS_GLEE *glee = (TS_GLEE*)ts->data; 737 PetscErrorCode ierr; 738 739 PetscFunctionBegin; 740 if (Ydot) { 741 if (dm && dm != ts->dm) { 742 ierr = DMGetNamedGlobalVector(dm,"TSGLEE_Ydot",Ydot);CHKERRQ(ierr); 743 } else *Ydot = glee->Ydot; 744 } 745 PetscFunctionReturn(0); 746 } 747 748 static PetscErrorCode TSGLEERestoreVecs(TS ts,DM dm,Vec *Ydot) 749 { 750 PetscErrorCode ierr; 751 752 PetscFunctionBegin; 753 if (Ydot) { 754 if (dm && dm != ts->dm) { 755 ierr = DMRestoreNamedGlobalVector(dm,"TSGLEE_Ydot",Ydot);CHKERRQ(ierr); 756 } 757 } 758 PetscFunctionReturn(0); 759 } 760 761 /* 762 This defines the nonlinear equation that is to be solved with SNES 763 */ 764 static PetscErrorCode SNESTSFormFunction_GLEE(SNES snes,Vec X,Vec F,TS ts) 765 { 766 TS_GLEE *glee = (TS_GLEE*)ts->data; 767 DM dm,dmsave; 768 Vec Ydot; 769 PetscReal shift = glee->scoeff / ts->time_step; 770 PetscErrorCode ierr; 771 772 PetscFunctionBegin; 773 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 774 ierr = TSGLEEGetVecs(ts,dm,&Ydot);CHKERRQ(ierr); 775 /* Set Ydot = shift*X */ 776 ierr = VecCopy(X,Ydot);CHKERRQ(ierr); 777 ierr = VecScale(Ydot,shift);CHKERRQ(ierr); 778 dmsave = ts->dm; 779 ts->dm = dm; 780 781 ierr = TSComputeIFunction(ts,glee->stage_time,X,Ydot,F,PETSC_FALSE);CHKERRQ(ierr); 782 783 ts->dm = dmsave; 784 ierr = TSGLEERestoreVecs(ts,dm,&Ydot);CHKERRQ(ierr); 785 PetscFunctionReturn(0); 786 } 787 788 static PetscErrorCode SNESTSFormJacobian_GLEE(SNES snes,Vec X,Mat A,Mat B,TS ts) 789 { 790 TS_GLEE *glee = (TS_GLEE*)ts->data; 791 DM dm,dmsave; 792 Vec Ydot; 793 PetscReal shift = glee->scoeff / ts->time_step; 794 PetscErrorCode ierr; 795 796 PetscFunctionBegin; 797 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 798 ierr = TSGLEEGetVecs(ts,dm,&Ydot);CHKERRQ(ierr); 799 /* glee->Ydot has already been computed in SNESTSFormFunction_GLEE (SNES guarantees this) */ 800 dmsave = ts->dm; 801 ts->dm = dm; 802 803 ierr = TSComputeIJacobian(ts,glee->stage_time,X,Ydot,shift,A,B,PETSC_FALSE);CHKERRQ(ierr); 804 805 ts->dm = dmsave; 806 ierr = TSGLEERestoreVecs(ts,dm,&Ydot);CHKERRQ(ierr); 807 PetscFunctionReturn(0); 808 } 809 810 static PetscErrorCode DMCoarsenHook_TSGLEE(DM fine,DM coarse,void *ctx) 811 { 812 PetscFunctionBegin; 813 PetscFunctionReturn(0); 814 } 815 816 static PetscErrorCode DMRestrictHook_TSGLEE(DM fine,Mat restrct,Vec rscale,Mat inject,DM coarse,void *ctx) 817 { 818 PetscFunctionBegin; 819 PetscFunctionReturn(0); 820 } 821 822 static PetscErrorCode DMSubDomainHook_TSGLEE(DM dm,DM subdm,void *ctx) 823 { 824 PetscFunctionBegin; 825 PetscFunctionReturn(0); 826 } 827 828 static PetscErrorCode DMSubDomainRestrictHook_TSGLEE(DM dm,VecScatter gscat,VecScatter lscat,DM subdm,void *ctx) 829 { 830 PetscFunctionBegin; 831 PetscFunctionReturn(0); 832 } 833 834 static PetscErrorCode TSSetUp_GLEE(TS ts) 835 { 836 TS_GLEE *glee = (TS_GLEE*)ts->data; 837 GLEETableau tab; 838 PetscInt s,r; 839 PetscErrorCode ierr; 840 DM dm; 841 842 PetscFunctionBegin; 843 if (!glee->tableau) { 844 ierr = TSGLEESetType(ts,TSGLEEDefaultType);CHKERRQ(ierr); 845 } 846 tab = glee->tableau; 847 s = tab->s; 848 r = tab->r; 849 ierr = VecDuplicateVecs(ts->vec_sol,r,&glee->Y);CHKERRQ(ierr); 850 ierr = VecDuplicateVecs(ts->vec_sol,r,&glee->X);CHKERRQ(ierr); 851 ierr = VecDuplicateVecs(ts->vec_sol,s,&glee->YStage);CHKERRQ(ierr); 852 ierr = VecDuplicateVecs(ts->vec_sol,s,&glee->YdotStage);CHKERRQ(ierr); 853 ierr = VecDuplicate(ts->vec_sol,&glee->Ydot);CHKERRQ(ierr); 854 ierr = VecDuplicate(ts->vec_sol,&glee->yGErr);CHKERRQ(ierr); 855 ierr = VecZeroEntries(glee->yGErr);CHKERRQ(ierr); 856 ierr = VecDuplicate(ts->vec_sol,&glee->W);CHKERRQ(ierr); 857 ierr = PetscMalloc2(s,&glee->swork,r,&glee->rwork);CHKERRQ(ierr); 858 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 859 ierr = DMCoarsenHookAdd(dm,DMCoarsenHook_TSGLEE,DMRestrictHook_TSGLEE,ts);CHKERRQ(ierr); 860 ierr = DMSubDomainHookAdd(dm,DMSubDomainHook_TSGLEE,DMSubDomainRestrictHook_TSGLEE,ts);CHKERRQ(ierr); 861 PetscFunctionReturn(0); 862 } 863 864 PetscErrorCode TSStartingMethod_GLEE(TS ts) 865 { 866 TS_GLEE *glee = (TS_GLEE*)ts->data; 867 GLEETableau tab = glee->tableau; 868 PetscInt r=tab->r,i; 869 PetscReal *S=tab->S; 870 PetscErrorCode ierr; 871 872 PetscFunctionBegin; 873 for (i=0; i<r; i++) { 874 ierr = VecZeroEntries(glee->Y[i]);CHKERRQ(ierr); 875 ierr = VecAXPY(glee->Y[i],S[i],ts->vec_sol);CHKERRQ(ierr); 876 } 877 878 PetscFunctionReturn(0); 879 } 880 881 /*------------------------------------------------------------*/ 882 883 static PetscErrorCode TSSetFromOptions_GLEE(PetscOptionItems *PetscOptionsObject,TS ts) 884 { 885 PetscErrorCode ierr; 886 char gleetype[256]; 887 888 PetscFunctionBegin; 889 ierr = PetscOptionsHead(PetscOptionsObject,"GLEE ODE solver options");CHKERRQ(ierr); 890 { 891 GLEETableauLink link; 892 PetscInt count,choice; 893 PetscBool flg; 894 const char **namelist; 895 896 ierr = PetscStrncpy(gleetype,TSGLEEDefaultType,sizeof(gleetype));CHKERRQ(ierr); 897 for (link=GLEETableauList,count=0; link; link=link->next,count++) ; 898 ierr = PetscMalloc1(count,(char***)&namelist);CHKERRQ(ierr); 899 for (link=GLEETableauList,count=0; link; link=link->next,count++) namelist[count] = link->tab.name; 900 ierr = PetscOptionsEList("-ts_glee_type","Family of GLEE method","TSGLEESetType",(const char*const*)namelist,count,gleetype,&choice,&flg);CHKERRQ(ierr); 901 ierr = TSGLEESetType(ts,flg ? namelist[choice] : gleetype);CHKERRQ(ierr); 902 ierr = PetscFree(namelist);CHKERRQ(ierr); 903 } 904 ierr = PetscOptionsTail();CHKERRQ(ierr); 905 PetscFunctionReturn(0); 906 } 907 908 static PetscErrorCode TSView_GLEE(TS ts,PetscViewer viewer) 909 { 910 TS_GLEE *glee = (TS_GLEE*)ts->data; 911 GLEETableau tab = glee->tableau; 912 PetscBool iascii; 913 PetscErrorCode ierr; 914 915 PetscFunctionBegin; 916 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 917 if (iascii) { 918 TSGLEEType gleetype; 919 char buf[512]; 920 ierr = TSGLEEGetType(ts,&gleetype);CHKERRQ(ierr); 921 ierr = PetscViewerASCIIPrintf(viewer," GLEE type %s\n",gleetype);CHKERRQ(ierr); 922 ierr = PetscFormatRealArray(buf,sizeof(buf),"% 8.6f",tab->s,tab->c);CHKERRQ(ierr); 923 ierr = PetscViewerASCIIPrintf(viewer," Abscissa c = %s\n",buf);CHKERRQ(ierr); 924 /* Note: print out r as well */ 925 } 926 PetscFunctionReturn(0); 927 } 928 929 static PetscErrorCode TSLoad_GLEE(TS ts,PetscViewer viewer) 930 { 931 PetscErrorCode ierr; 932 SNES snes; 933 TSAdapt tsadapt; 934 935 PetscFunctionBegin; 936 ierr = TSGetAdapt(ts,&tsadapt);CHKERRQ(ierr); 937 ierr = TSAdaptLoad(tsadapt,viewer);CHKERRQ(ierr); 938 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 939 ierr = SNESLoad(snes,viewer);CHKERRQ(ierr); 940 /* function and Jacobian context for SNES when used with TS is always ts object */ 941 ierr = SNESSetFunction(snes,NULL,NULL,ts);CHKERRQ(ierr); 942 ierr = SNESSetJacobian(snes,NULL,NULL,NULL,ts);CHKERRQ(ierr); 943 PetscFunctionReturn(0); 944 } 945 946 /*@C 947 TSGLEESetType - Set the type of GLEE scheme 948 949 Logically collective 950 951 Input Parameters: 952 + ts - timestepping context 953 - gleetype - type of GLEE-scheme 954 955 Level: intermediate 956 957 .seealso: TSGLEEGetType(), TSGLEE 958 @*/ 959 PetscErrorCode TSGLEESetType(TS ts,TSGLEEType gleetype) 960 { 961 PetscErrorCode ierr; 962 963 PetscFunctionBegin; 964 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 965 PetscValidCharPointer(gleetype,2); 966 ierr = PetscTryMethod(ts,"TSGLEESetType_C",(TS,TSGLEEType),(ts,gleetype));CHKERRQ(ierr); 967 PetscFunctionReturn(0); 968 } 969 970 /*@C 971 TSGLEEGetType - Get the type of GLEE scheme 972 973 Logically collective 974 975 Input Parameter: 976 . ts - timestepping context 977 978 Output Parameter: 979 . gleetype - type of GLEE-scheme 980 981 Level: intermediate 982 983 .seealso: TSGLEESetType() 984 @*/ 985 PetscErrorCode TSGLEEGetType(TS ts,TSGLEEType *gleetype) 986 { 987 PetscErrorCode ierr; 988 989 PetscFunctionBegin; 990 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 991 ierr = PetscUseMethod(ts,"TSGLEEGetType_C",(TS,TSGLEEType*),(ts,gleetype));CHKERRQ(ierr); 992 PetscFunctionReturn(0); 993 } 994 995 PetscErrorCode TSGLEEGetType_GLEE(TS ts,TSGLEEType *gleetype) 996 { 997 TS_GLEE *glee = (TS_GLEE*)ts->data; 998 PetscErrorCode ierr; 999 1000 PetscFunctionBegin; 1001 if (!glee->tableau) { 1002 ierr = TSGLEESetType(ts,TSGLEEDefaultType);CHKERRQ(ierr); 1003 } 1004 *gleetype = glee->tableau->name; 1005 PetscFunctionReturn(0); 1006 } 1007 PetscErrorCode TSGLEESetType_GLEE(TS ts,TSGLEEType gleetype) 1008 { 1009 TS_GLEE *glee = (TS_GLEE*)ts->data; 1010 PetscErrorCode ierr; 1011 PetscBool match; 1012 GLEETableauLink link; 1013 1014 PetscFunctionBegin; 1015 if (glee->tableau) { 1016 ierr = PetscStrcmp(glee->tableau->name,gleetype,&match);CHKERRQ(ierr); 1017 if (match) PetscFunctionReturn(0); 1018 } 1019 for (link = GLEETableauList; link; link=link->next) { 1020 ierr = PetscStrcmp(link->tab.name,gleetype,&match);CHKERRQ(ierr); 1021 if (match) { 1022 ierr = TSReset_GLEE(ts);CHKERRQ(ierr); 1023 glee->tableau = &link->tab; 1024 PetscFunctionReturn(0); 1025 } 1026 } 1027 SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_UNKNOWN_TYPE,"Could not find '%s'",gleetype); 1028 } 1029 1030 static PetscErrorCode TSGetStages_GLEE(TS ts,PetscInt *ns,Vec **Y) 1031 { 1032 TS_GLEE *glee = (TS_GLEE*)ts->data; 1033 1034 PetscFunctionBegin; 1035 if (ns) *ns = glee->tableau->s; 1036 if (Y) *Y = glee->YStage; 1037 PetscFunctionReturn(0); 1038 } 1039 1040 PetscErrorCode TSGetSolutionComponents_GLEE(TS ts,PetscInt *n,Vec *Y) 1041 { 1042 TS_GLEE *glee = (TS_GLEE*)ts->data; 1043 GLEETableau tab = glee->tableau; 1044 PetscErrorCode ierr; 1045 1046 PetscFunctionBegin; 1047 if (!Y) *n = tab->r; 1048 else { 1049 if ((*n >= 0) && (*n < tab->r)) { 1050 ierr = VecCopy(glee->Y[*n],*Y);CHKERRQ(ierr); 1051 } else SETERRQ3(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_OUTOFRANGE,"Second argument (%d) out of range[%d,%d].",*n,0,tab->r-1); 1052 } 1053 PetscFunctionReturn(0); 1054 } 1055 1056 PetscErrorCode TSGetAuxSolution_GLEE(TS ts,Vec *X) 1057 { 1058 TS_GLEE *glee = (TS_GLEE*)ts->data; 1059 GLEETableau tab = glee->tableau; 1060 PetscReal *F = tab->Fembed; 1061 PetscInt r = tab->r; 1062 Vec *Y = glee->Y; 1063 PetscScalar *wr = glee->rwork; 1064 PetscInt i; 1065 PetscErrorCode ierr; 1066 1067 PetscFunctionBegin; 1068 ierr = VecZeroEntries(*X);CHKERRQ(ierr); 1069 for (i=0; i<r; i++) wr[i] = F[i]; 1070 ierr = VecMAXPY((*X),r,wr,Y);CHKERRQ(ierr); 1071 PetscFunctionReturn(0); 1072 } 1073 1074 PetscErrorCode TSGetTimeError_GLEE(TS ts,PetscInt n,Vec *X) 1075 { 1076 TS_GLEE *glee = (TS_GLEE*)ts->data; 1077 GLEETableau tab = glee->tableau; 1078 PetscReal *F = tab->Ferror; 1079 PetscInt r = tab->r; 1080 Vec *Y = glee->Y; 1081 PetscScalar *wr = glee->rwork; 1082 PetscInt i; 1083 PetscErrorCode ierr; 1084 1085 PetscFunctionBegin; 1086 ierr = VecZeroEntries(*X);CHKERRQ(ierr); 1087 if (n==0) { 1088 for (i=0; i<r; i++) wr[i] = F[i]; 1089 ierr = VecMAXPY((*X),r,wr,Y);CHKERRQ(ierr); 1090 } else if (n==-1) { 1091 *X=glee->yGErr; 1092 } 1093 PetscFunctionReturn(0); 1094 } 1095 1096 PetscErrorCode TSSetTimeError_GLEE(TS ts,Vec X) 1097 { 1098 TS_GLEE *glee = (TS_GLEE*)ts->data; 1099 GLEETableau tab = glee->tableau; 1100 PetscReal *S = tab->Serror; 1101 PetscInt r = tab->r,i; 1102 Vec *Y = glee->Y; 1103 PetscErrorCode ierr; 1104 1105 PetscFunctionBegin; 1106 if (r != 2) SETERRQ2(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"TSSetTimeError_GLEE not supported for '%s' with r=%D.",tab->name,tab->r); 1107 for (i=1; i<r; i++) { 1108 ierr = VecCopy(ts->vec_sol,Y[i]);CHKERRQ(ierr); 1109 ierr = VecAXPBY(Y[i],S[0],S[1],X);CHKERRQ(ierr); 1110 ierr = VecCopy(X,glee->yGErr);CHKERRQ(ierr); 1111 } 1112 PetscFunctionReturn(0); 1113 } 1114 1115 static PetscErrorCode TSDestroy_GLEE(TS ts) 1116 { 1117 PetscErrorCode ierr; 1118 1119 PetscFunctionBegin; 1120 ierr = TSReset_GLEE(ts);CHKERRQ(ierr); 1121 if (ts->dm) { 1122 ierr = DMCoarsenHookRemove(ts->dm,DMCoarsenHook_TSGLEE,DMRestrictHook_TSGLEE,ts);CHKERRQ(ierr); 1123 ierr = DMSubDomainHookRemove(ts->dm,DMSubDomainHook_TSGLEE,DMSubDomainRestrictHook_TSGLEE,ts);CHKERRQ(ierr); 1124 } 1125 ierr = PetscFree(ts->data);CHKERRQ(ierr); 1126 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSGLEEGetType_C",NULL);CHKERRQ(ierr); 1127 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSGLEESetType_C",NULL);CHKERRQ(ierr); 1128 PetscFunctionReturn(0); 1129 } 1130 1131 /* ------------------------------------------------------------ */ 1132 /*MC 1133 TSGLEE - ODE and DAE solver using General Linear with Error Estimation schemes 1134 1135 The user should provide the right hand side of the equation 1136 using TSSetRHSFunction(). 1137 1138 Notes: 1139 The default is TSGLEE35, it can be changed with TSGLEESetType() or -ts_glee_type 1140 1141 Level: beginner 1142 1143 .seealso: TSCreate(), TS, TSSetType(), TSGLEESetType(), TSGLEEGetType(), 1144 TSGLEE23, TTSGLEE24, TSGLEE35, TSGLEE25I, TSGLEEEXRK2A, 1145 TSGLEERK32G1, TSGLEERK285EX, TSGLEEType, TSGLEERegister() 1146 1147 M*/ 1148 PETSC_EXTERN PetscErrorCode TSCreate_GLEE(TS ts) 1149 { 1150 TS_GLEE *th; 1151 PetscErrorCode ierr; 1152 1153 PetscFunctionBegin; 1154 ierr = TSGLEEInitializePackage();CHKERRQ(ierr); 1155 1156 ts->ops->reset = TSReset_GLEE; 1157 ts->ops->destroy = TSDestroy_GLEE; 1158 ts->ops->view = TSView_GLEE; 1159 ts->ops->load = TSLoad_GLEE; 1160 ts->ops->setup = TSSetUp_GLEE; 1161 ts->ops->step = TSStep_GLEE; 1162 ts->ops->interpolate = TSInterpolate_GLEE; 1163 ts->ops->evaluatestep = TSEvaluateStep_GLEE; 1164 ts->ops->setfromoptions = TSSetFromOptions_GLEE; 1165 ts->ops->getstages = TSGetStages_GLEE; 1166 ts->ops->snesfunction = SNESTSFormFunction_GLEE; 1167 ts->ops->snesjacobian = SNESTSFormJacobian_GLEE; 1168 ts->ops->getsolutioncomponents = TSGetSolutionComponents_GLEE; 1169 ts->ops->getauxsolution = TSGetAuxSolution_GLEE; 1170 ts->ops->gettimeerror = TSGetTimeError_GLEE; 1171 ts->ops->settimeerror = TSSetTimeError_GLEE; 1172 ts->ops->startingmethod = TSStartingMethod_GLEE; 1173 ts->default_adapt_type = TSADAPTGLEE; 1174 1175 ts->usessnes = PETSC_TRUE; 1176 1177 ierr = PetscNewLog(ts,&th);CHKERRQ(ierr); 1178 ts->data = (void*)th; 1179 1180 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSGLEEGetType_C",TSGLEEGetType_GLEE);CHKERRQ(ierr); 1181 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSGLEESetType_C",TSGLEESetType_GLEE);CHKERRQ(ierr); 1182 PetscFunctionReturn(0); 1183 } 1184