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