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 = PetscMalloc(sizeof(*link),&link);CHKERRQ(ierr); 477 ierr = PetscMemzero(link,sizeof(*link));CHKERRQ(ierr); 478 t = &link->tab; 479 ierr = PetscStrallocpy(name,&t->name);CHKERRQ(ierr); 480 t->order = order; 481 t->s = s; 482 t->r = r; 483 t->gamma = gamma; 484 ierr = PetscMalloc5(s*s,&t->A,r*r,&t->V,s,&t->c,r*s,&t->B,s*r,&t->U); CHKERRQ(ierr); 485 ierr = PetscMalloc2(r,&t->S,r,&t->F); CHKERRQ(ierr); 486 ierr = PetscMemcpy(t->A,A,s*s*sizeof(A[0]));CHKERRQ(ierr); 487 ierr = PetscMemcpy(t->B,B,r*s*sizeof(B[0]));CHKERRQ(ierr); 488 ierr = PetscMemcpy(t->U,U,s*r*sizeof(U[0]));CHKERRQ(ierr); 489 ierr = PetscMemcpy(t->V,V,r*r*sizeof(V[0]));CHKERRQ(ierr); 490 ierr = PetscMemcpy(t->S,S,r *sizeof(S[0]));CHKERRQ(ierr); 491 ierr = PetscMemcpy(t->F,F,r *sizeof(F[0]));CHKERRQ(ierr); 492 if (c) { 493 ierr = PetscMemcpy(t->c,c,s*sizeof(c[0])); CHKERRQ(ierr); 494 } else { 495 for (i=0; i<s; i++) for (j=0,t->c[i]=0; j<s; j++) t->c[i] += A[i*s+j]; 496 } 497 ierr = PetscMalloc1(r,&t->Fembed);CHKERRQ(ierr); 498 ierr = PetscMalloc1(r,&t->Ferror);CHKERRQ(ierr); 499 ierr = PetscMalloc1(r,&t->Serror);CHKERRQ(ierr); 500 ierr = PetscMemcpy(t->Fembed,Fembed,r*sizeof(Fembed[0]));CHKERRQ(ierr); 501 ierr = PetscMemcpy(t->Ferror,Ferror,r*sizeof(Ferror[0]));CHKERRQ(ierr); 502 ierr = PetscMemcpy(t->Serror,Serror,r*sizeof(Serror[0]));CHKERRQ(ierr); 503 t->pinterp = pinterp; 504 ierr = PetscMalloc1(s*pinterp,&t->binterp);CHKERRQ(ierr); 505 ierr = PetscMemcpy(t->binterp,binterp,s*pinterp*sizeof(binterp[0]));CHKERRQ(ierr); 506 507 link->next = GLEETableauList; 508 GLEETableauList = link; 509 PetscFunctionReturn(0); 510 } 511 512 static PetscErrorCode TSEvaluateStep_GLEE(TS ts,PetscInt order,Vec X,PetscBool *done) 513 { 514 TS_GLEE *glee = (TS_GLEE*) ts->data; 515 GLEETableau tab = glee->tableau; 516 PetscReal h, *B = tab->B, *V = tab->V, 517 *F = tab->F, 518 *Fembed = tab->Fembed; 519 PetscInt s = tab->s, r = tab->r, i, j; 520 Vec *Y = glee->Y, *YdotStage = glee->YdotStage; 521 PetscScalar *ws = glee->swork, *wr = glee->rwork; 522 PetscErrorCode ierr; 523 524 PetscFunctionBegin; 525 526 switch (glee->status) { 527 case TS_STEP_INCOMPLETE: 528 case TS_STEP_PENDING: 529 h = ts->time_step; break; 530 case TS_STEP_COMPLETE: 531 h = ts->ptime - ts->ptime_prev; break; 532 default: SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_PLIB,"Invalid TSStepStatus"); 533 } 534 535 if (order == tab->order) { 536 537 /* Note: Irrespective of whether status is TS_STEP_INCOMPLETE 538 or TS_STEP_COMPLETE, glee->X has the solution at the 539 beginning of the time step. So no need to roll-back. 540 */ 541 if (glee->status == TS_STEP_INCOMPLETE) { 542 for (i=0; i<r; i++) { 543 ierr = VecZeroEntries(Y[i]); CHKERRQ(ierr); 544 for (j=0; j<r; j++) wr[j] = V[i*r+j]; 545 ierr = VecMAXPY(Y[i],r,wr,glee->X); CHKERRQ(ierr); 546 for (j=0; j<s; j++) ws[j] = h*B[i*s+j]; 547 ierr = VecMAXPY(Y[i],s,ws,YdotStage); CHKERRQ(ierr); 548 } 549 ierr = VecZeroEntries(X);CHKERRQ(ierr); 550 for (j=0; j<r; j++) wr[j] = F[j]; 551 ierr = VecMAXPY(X,r,wr,Y);CHKERRQ(ierr); 552 } else {ierr = VecCopy(ts->vec_sol,X);CHKERRQ(ierr);} 553 PetscFunctionReturn(0); 554 555 } else if (order == tab->order-1) { 556 557 /* Complete with the embedded method (Fembed) */ 558 for (i=0; i<r; i++) { 559 ierr = VecZeroEntries(Y[i]); CHKERRQ(ierr); 560 for (j=0; j<r; j++) wr[j] = V[i*r+j]; 561 ierr = VecMAXPY(Y[i],r,wr,glee->X); CHKERRQ(ierr); 562 for (j=0; j<s; j++) ws[j] = h*B[i*s+j]; 563 ierr = VecMAXPY(Y[i],s,ws,YdotStage); CHKERRQ(ierr); 564 } 565 ierr = VecZeroEntries(X);CHKERRQ(ierr); 566 for (j=0; j<r; j++) wr[j] = Fembed[j]; 567 ierr = VecMAXPY(X,r,wr,Y);CHKERRQ(ierr); 568 569 if (done) *done = PETSC_TRUE; 570 PetscFunctionReturn(0); 571 } 572 if (done) *done = PETSC_FALSE; 573 else SETERRQ3(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"GLEE '%s' of order %D cannot evaluate step at order %D",tab->name,tab->order,order); 574 PetscFunctionReturn(0); 575 } 576 577 static PetscErrorCode TSStep_GLEE(TS ts) 578 { 579 TS_GLEE *glee = (TS_GLEE*)ts->data; 580 GLEETableau tab = glee->tableau; 581 const PetscInt s = tab->s, r = tab->r; 582 PetscReal *A = tab->A, *U = tab->U, 583 *F = tab->F, 584 *c = tab->c; 585 Vec *Y = glee->Y, *X = glee->X, 586 *YStage = glee->YStage, 587 *YdotStage = glee->YdotStage, 588 W = glee->W; 589 SNES snes; 590 PetscScalar *ws = glee->swork, *wr = glee->rwork; 591 TSAdapt adapt; 592 PetscInt i,j,reject,next_scheme,its,lits; 593 PetscReal next_time_step; 594 PetscReal t; 595 PetscBool accept; 596 PetscErrorCode ierr; 597 598 PetscFunctionBegin; 599 ierr = PetscCitationsRegister(citation,&cited);CHKERRQ(ierr); 600 601 for (i=0; i<r; i++) { ierr = VecCopy(Y[i],X[i]); CHKERRQ(ierr); } 602 603 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 604 next_time_step = ts->time_step; 605 t = ts->ptime; 606 accept = PETSC_TRUE; 607 glee->status = TS_STEP_INCOMPLETE; 608 609 for (reject=0; reject<ts->max_reject && !ts->reason; reject++,ts->reject++) { 610 611 PetscReal h = ts->time_step; 612 ierr = TSPreStep(ts);CHKERRQ(ierr); 613 614 for (i=0; i<s; i++) { 615 616 glee->stage_time = t + h*c[i]; 617 ierr = TSPreStage(ts,glee->stage_time); CHKERRQ(ierr); 618 619 if (A[i*s+i] == 0) { /* Explicit stage */ 620 ierr = VecZeroEntries(YStage[i]);CHKERRQ(ierr); 621 for (j=0; j<r; j++) wr[j] = U[i*r+j]; 622 ierr = VecMAXPY(YStage[i],r,wr,X);CHKERRQ(ierr); 623 for (j=0; j<i; j++) ws[j] = h*A[i*s+j]; 624 ierr = VecMAXPY(YStage[i],i,ws,YdotStage);CHKERRQ(ierr); 625 } else { /* Implicit stage */ 626 glee->scoeff = 1.0/A[i*s+i]; 627 /* compute right-hand-side */ 628 ierr = VecZeroEntries(W);CHKERRQ(ierr); 629 for (j=0; j<r; j++) wr[j] = U[i*r+j]; 630 ierr = VecMAXPY(W,r,wr,X);CHKERRQ(ierr); 631 for (j=0; j<i; j++) ws[j] = h*A[i*s+j]; 632 ierr = VecMAXPY(W,i,ws,YdotStage);CHKERRQ(ierr); 633 ierr = VecScale(W,glee->scoeff/h);CHKERRQ(ierr); 634 /* set initial guess */ 635 ierr = VecCopy(i>0 ? YStage[i-1] : ts->vec_sol,YStage[i]);CHKERRQ(ierr); 636 /* solve for this stage */ 637 ierr = SNESSolve(snes,W,YStage[i]);CHKERRQ(ierr); 638 ierr = SNESGetIterationNumber(snes,&its);CHKERRQ(ierr); 639 ierr = SNESGetLinearSolveIterations(snes,&lits);CHKERRQ(ierr); 640 ts->snes_its += its; ts->ksp_its += lits; 641 } 642 ierr = TSGetAdapt(ts,&adapt);CHKERRQ(ierr); 643 ierr = TSAdaptCheckStage(adapt,ts,glee->stage_time,YStage[i],&accept);CHKERRQ(ierr); 644 if (!accept) goto reject_step; 645 ierr = TSPostStage(ts,glee->stage_time,i,YStage); CHKERRQ(ierr); 646 ierr = TSComputeRHSFunction(ts,t+h*c[i],YStage[i],YdotStage[i]);CHKERRQ(ierr); 647 } 648 ierr = TSEvaluateStep(ts,tab->order,ts->vec_sol,NULL);CHKERRQ(ierr); 649 glee->status = TS_STEP_PENDING; 650 651 /* Register only the current method as a candidate because we're not supporting multiple candidates yet. */ 652 ierr = TSGetAdapt(ts,&adapt);CHKERRQ(ierr); 653 ierr = TSAdaptCandidatesClear(adapt);CHKERRQ(ierr); 654 ierr = TSAdaptCandidateAdd(adapt,tab->name,tab->order,1,tab->ccfl,(PetscReal)tab->s,PETSC_TRUE);CHKERRQ(ierr); 655 ierr = TSAdaptChoose(adapt,ts,ts->time_step,&next_scheme,&next_time_step,&accept);CHKERRQ(ierr); 656 if (accept) { 657 /* ignore next_scheme for now */ 658 ts->ptime += ts->time_step; 659 ts->time_step = next_time_step; 660 glee->status = TS_STEP_COMPLETE; 661 /* compute and store the global error */ 662 /* Note: this is not needed if TSAdaptGLEE is not used */ 663 ierr = TSGetTimeError(ts,0,&(glee->yGErr));CHKERRQ(ierr); 664 ierr = PetscObjectComposedDataSetReal((PetscObject)ts->vec_sol,explicit_stage_time_id,ts->ptime);CHKERRQ(ierr); 665 break; 666 } else { /* Roll back the current step */ 667 for (j=0; j<r; j++) wr[j] = F[j]; 668 ierr = VecMAXPY(ts->vec_sol,r,wr,X);CHKERRQ(ierr); 669 ts->time_step = next_time_step; 670 glee->status = TS_STEP_INCOMPLETE; 671 } 672 reject_step: continue; 673 } 674 if (glee->status != TS_STEP_COMPLETE && !ts->reason) ts->reason = TS_DIVERGED_STEP_REJECTED; 675 PetscFunctionReturn(0); 676 } 677 678 static PetscErrorCode TSInterpolate_GLEE(TS ts,PetscReal itime,Vec X) 679 { 680 TS_GLEE *glee = (TS_GLEE*)ts->data; 681 PetscInt s=glee->tableau->s, pinterp=glee->tableau->pinterp,i,j; 682 PetscReal h,tt,t; 683 PetscScalar *b; 684 const PetscReal *B = glee->tableau->binterp; 685 PetscErrorCode ierr; 686 687 PetscFunctionBegin; 688 if (!B) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"TSGLEE %s does not have an interpolation formula",glee->tableau->name); 689 switch (glee->status) { 690 case TS_STEP_INCOMPLETE: 691 case TS_STEP_PENDING: 692 h = ts->time_step; 693 t = (itime - ts->ptime)/h; 694 break; 695 case TS_STEP_COMPLETE: 696 h = ts->ptime - ts->ptime_prev; 697 t = (itime - ts->ptime)/h + 1; /* In the interval [0,1] */ 698 break; 699 default: SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_PLIB,"Invalid TSStepStatus"); 700 } 701 ierr = PetscMalloc1(s,&b);CHKERRQ(ierr); 702 for (i=0; i<s; i++) b[i] = 0; 703 for (j=0,tt=t; j<pinterp; j++,tt*=t) { 704 for (i=0; i<s; i++) { 705 b[i] += h * B[i*pinterp+j] * tt; 706 } 707 } 708 ierr = VecCopy(glee->YStage[0],X);CHKERRQ(ierr); 709 ierr = VecMAXPY(X,s,b,glee->YdotStage);CHKERRQ(ierr); 710 ierr = PetscFree(b);CHKERRQ(ierr); 711 PetscFunctionReturn(0); 712 } 713 714 /*------------------------------------------------------------*/ 715 static PetscErrorCode TSReset_GLEE(TS ts) 716 { 717 TS_GLEE *glee = (TS_GLEE*)ts->data; 718 PetscInt s, r; 719 PetscErrorCode ierr; 720 721 PetscFunctionBegin; 722 if (!glee->tableau) PetscFunctionReturn(0); 723 s = glee->tableau->s; 724 r = glee->tableau->r; 725 ierr = VecDestroyVecs(r,&glee->Y);CHKERRQ(ierr); 726 ierr = VecDestroyVecs(r,&glee->X);CHKERRQ(ierr); 727 ierr = VecDestroyVecs(s,&glee->YStage);CHKERRQ(ierr); 728 ierr = VecDestroyVecs(s,&glee->YdotStage);CHKERRQ(ierr); 729 ierr = VecDestroy(&glee->Ydot);CHKERRQ(ierr); 730 ierr = VecDestroy(&glee->yGErr);CHKERRQ(ierr); 731 ierr = VecDestroy(&glee->W);CHKERRQ(ierr); 732 ierr = PetscFree2(glee->swork,glee->rwork);CHKERRQ(ierr); 733 PetscFunctionReturn(0); 734 } 735 736 static PetscErrorCode TSGLEEGetVecs(TS ts,DM dm,Vec *Ydot) 737 { 738 TS_GLEE *glee = (TS_GLEE*)ts->data; 739 PetscErrorCode ierr; 740 741 PetscFunctionBegin; 742 if (Ydot) { 743 if (dm && dm != ts->dm) { 744 ierr = DMGetNamedGlobalVector(dm,"TSGLEE_Ydot",Ydot);CHKERRQ(ierr); 745 } else *Ydot = glee->Ydot; 746 } 747 PetscFunctionReturn(0); 748 } 749 750 751 static PetscErrorCode TSGLEERestoreVecs(TS ts,DM dm,Vec *Ydot) 752 { 753 PetscErrorCode ierr; 754 755 PetscFunctionBegin; 756 if (Ydot) { 757 if (dm && dm != ts->dm) { 758 ierr = DMRestoreNamedGlobalVector(dm,"TSGLEE_Ydot",Ydot);CHKERRQ(ierr); 759 } 760 } 761 PetscFunctionReturn(0); 762 } 763 764 /* 765 This defines the nonlinear equation that is to be solved with SNES 766 */ 767 static PetscErrorCode SNESTSFormFunction_GLEE(SNES snes,Vec X,Vec F,TS ts) 768 { 769 TS_GLEE *glee = (TS_GLEE*)ts->data; 770 DM dm,dmsave; 771 Vec Ydot; 772 PetscReal shift = glee->scoeff / ts->time_step; 773 PetscErrorCode ierr; 774 775 PetscFunctionBegin; 776 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 777 ierr = TSGLEEGetVecs(ts,dm,&Ydot);CHKERRQ(ierr); 778 /* Set Ydot = shift*X */ 779 ierr = VecCopy(X,Ydot);CHKERRQ(ierr); 780 ierr = VecScale(Ydot,shift);CHKERRQ(ierr); 781 dmsave = ts->dm; 782 ts->dm = dm; 783 784 ierr = TSComputeIFunction(ts,glee->stage_time,X,Ydot,F,PETSC_FALSE);CHKERRQ(ierr); 785 786 ts->dm = dmsave; 787 ierr = TSGLEERestoreVecs(ts,dm,&Ydot);CHKERRQ(ierr); 788 PetscFunctionReturn(0); 789 } 790 791 static PetscErrorCode SNESTSFormJacobian_GLEE(SNES snes,Vec X,Mat A,Mat B,TS ts) 792 { 793 TS_GLEE *glee = (TS_GLEE*)ts->data; 794 DM dm,dmsave; 795 Vec Ydot; 796 PetscReal shift = glee->scoeff / ts->time_step; 797 PetscErrorCode ierr; 798 799 PetscFunctionBegin; 800 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 801 ierr = TSGLEEGetVecs(ts,dm,&Ydot);CHKERRQ(ierr); 802 /* glee->Ydot has already been computed in SNESTSFormFunction_GLEE (SNES guarantees this) */ 803 dmsave = ts->dm; 804 ts->dm = dm; 805 806 ierr = TSComputeIJacobian(ts,glee->stage_time,X,Ydot,shift,A,B,PETSC_FALSE);CHKERRQ(ierr); 807 808 ts->dm = dmsave; 809 ierr = TSGLEERestoreVecs(ts,dm,&Ydot);CHKERRQ(ierr); 810 PetscFunctionReturn(0); 811 } 812 813 static PetscErrorCode DMCoarsenHook_TSGLEE(DM fine,DM coarse,void *ctx) 814 { 815 PetscFunctionBegin; 816 PetscFunctionReturn(0); 817 } 818 819 static PetscErrorCode DMRestrictHook_TSGLEE(DM fine,Mat restrct,Vec rscale,Mat inject,DM coarse,void *ctx) 820 { 821 PetscFunctionBegin; 822 PetscFunctionReturn(0); 823 } 824 825 826 static PetscErrorCode DMSubDomainHook_TSGLEE(DM dm,DM subdm,void *ctx) 827 { 828 PetscFunctionBegin; 829 PetscFunctionReturn(0); 830 } 831 832 static PetscErrorCode DMSubDomainRestrictHook_TSGLEE(DM dm,VecScatter gscat,VecScatter lscat,DM subdm,void *ctx) 833 { 834 PetscFunctionBegin; 835 PetscFunctionReturn(0); 836 } 837 838 static PetscErrorCode TSSetUp_GLEE(TS ts) 839 { 840 TS_GLEE *glee = (TS_GLEE*)ts->data; 841 GLEETableau tab; 842 PetscInt s,r; 843 PetscErrorCode ierr; 844 DM dm; 845 846 PetscFunctionBegin; 847 if (!glee->tableau) { 848 ierr = TSGLEESetType(ts,TSGLEEDefaultType);CHKERRQ(ierr); 849 } 850 tab = glee->tableau; 851 s = tab->s; 852 r = tab->r; 853 ierr = VecDuplicateVecs(ts->vec_sol,r,&glee->Y);CHKERRQ(ierr); 854 ierr = VecDuplicateVecs(ts->vec_sol,r,&glee->X);CHKERRQ(ierr); 855 ierr = VecDuplicateVecs(ts->vec_sol,s,&glee->YStage);CHKERRQ(ierr); 856 ierr = VecDuplicateVecs(ts->vec_sol,s,&glee->YdotStage);CHKERRQ(ierr); 857 ierr = VecDuplicate(ts->vec_sol,&glee->Ydot);CHKERRQ(ierr); 858 ierr = VecDuplicate(ts->vec_sol,&glee->yGErr);CHKERRQ(ierr); 859 ierr = VecZeroEntries(glee->yGErr);CHKERRQ(ierr); 860 ierr = VecDuplicate(ts->vec_sol,&glee->W);CHKERRQ(ierr); 861 ierr = PetscMalloc2(s,&glee->swork,r,&glee->rwork); CHKERRQ(ierr); 862 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 863 ierr = DMCoarsenHookAdd(dm,DMCoarsenHook_TSGLEE,DMRestrictHook_TSGLEE,ts);CHKERRQ(ierr); 864 ierr = DMSubDomainHookAdd(dm,DMSubDomainHook_TSGLEE,DMSubDomainRestrictHook_TSGLEE,ts);CHKERRQ(ierr); 865 PetscFunctionReturn(0); 866 } 867 868 PetscErrorCode TSStartingMethod_GLEE(TS ts) 869 { 870 TS_GLEE *glee = (TS_GLEE*)ts->data; 871 GLEETableau tab = glee->tableau; 872 PetscInt r=tab->r,i; 873 PetscReal *S=tab->S; 874 PetscErrorCode ierr; 875 876 PetscFunctionBegin; 877 for (i=0; i<r; i++) { 878 ierr = VecZeroEntries(glee->Y[i]);CHKERRQ(ierr); 879 ierr = VecAXPY(glee->Y[i],S[i],ts->vec_sol);CHKERRQ(ierr); 880 } 881 882 PetscFunctionReturn(0); 883 } 884 885 886 /*------------------------------------------------------------*/ 887 888 static PetscErrorCode TSSetFromOptions_GLEE(PetscOptionItems *PetscOptionsObject,TS ts) 889 { 890 PetscErrorCode ierr; 891 char gleetype[256]; 892 893 PetscFunctionBegin; 894 ierr = PetscOptionsHead(PetscOptionsObject,"GLEE ODE solver options");CHKERRQ(ierr); 895 { 896 GLEETableauLink link; 897 PetscInt count,choice; 898 PetscBool flg; 899 const char **namelist; 900 901 ierr = PetscStrncpy(gleetype,TSGLEEDefaultType,sizeof(gleetype));CHKERRQ(ierr); 902 for (link=GLEETableauList,count=0; link; link=link->next,count++) ; 903 ierr = PetscMalloc1(count,(char***)&namelist);CHKERRQ(ierr); 904 for (link=GLEETableauList,count=0; link; link=link->next,count++) namelist[count] = link->tab.name; 905 ierr = PetscOptionsEList("-ts_glee_type","Family of GLEE method","TSGLEESetType",(const char*const*)namelist,count,gleetype,&choice,&flg);CHKERRQ(ierr); 906 ierr = TSGLEESetType(ts,flg ? namelist[choice] : gleetype);CHKERRQ(ierr); 907 ierr = PetscFree(namelist);CHKERRQ(ierr); 908 } 909 ierr = PetscOptionsTail();CHKERRQ(ierr); 910 PetscFunctionReturn(0); 911 } 912 913 static PetscErrorCode TSView_GLEE(TS ts,PetscViewer viewer) 914 { 915 TS_GLEE *glee = (TS_GLEE*)ts->data; 916 GLEETableau tab = glee->tableau; 917 PetscBool iascii; 918 PetscErrorCode ierr; 919 920 PetscFunctionBegin; 921 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 922 if (iascii) { 923 TSGLEEType gleetype; 924 char buf[512]; 925 ierr = TSGLEEGetType(ts,&gleetype);CHKERRQ(ierr); 926 ierr = PetscViewerASCIIPrintf(viewer," GLEE type %s\n",gleetype);CHKERRQ(ierr); 927 ierr = PetscFormatRealArray(buf,sizeof(buf),"% 8.6f",tab->s,tab->c);CHKERRQ(ierr); 928 ierr = PetscViewerASCIIPrintf(viewer," Abscissa c = %s\n",buf);CHKERRQ(ierr); 929 /* Note: print out r as well */ 930 } 931 PetscFunctionReturn(0); 932 } 933 934 static PetscErrorCode TSLoad_GLEE(TS ts,PetscViewer viewer) 935 { 936 PetscErrorCode ierr; 937 SNES snes; 938 TSAdapt tsadapt; 939 940 PetscFunctionBegin; 941 ierr = TSGetAdapt(ts,&tsadapt);CHKERRQ(ierr); 942 ierr = TSAdaptLoad(tsadapt,viewer);CHKERRQ(ierr); 943 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 944 ierr = SNESLoad(snes,viewer);CHKERRQ(ierr); 945 /* function and Jacobian context for SNES when used with TS is always ts object */ 946 ierr = SNESSetFunction(snes,NULL,NULL,ts);CHKERRQ(ierr); 947 ierr = SNESSetJacobian(snes,NULL,NULL,NULL,ts);CHKERRQ(ierr); 948 PetscFunctionReturn(0); 949 } 950 951 /*@C 952 TSGLEESetType - Set the type of GLEE scheme 953 954 Logically collective 955 956 Input Parameter: 957 + ts - timestepping context 958 - gleetype - type of GLEE-scheme 959 960 Level: intermediate 961 962 .seealso: TSGLEEGetType(), TSGLEE 963 @*/ 964 PetscErrorCode TSGLEESetType(TS ts,TSGLEEType gleetype) 965 { 966 PetscErrorCode ierr; 967 968 PetscFunctionBegin; 969 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 970 PetscValidCharPointer(gleetype,2); 971 ierr = PetscTryMethod(ts,"TSGLEESetType_C",(TS,TSGLEEType),(ts,gleetype));CHKERRQ(ierr); 972 PetscFunctionReturn(0); 973 } 974 975 /*@C 976 TSGLEEGetType - Get the type of GLEE scheme 977 978 Logically collective 979 980 Input Parameter: 981 . ts - timestepping context 982 983 Output Parameter: 984 . gleetype - type of GLEE-scheme 985 986 Level: intermediate 987 988 .seealso: TSGLEESetType() 989 @*/ 990 PetscErrorCode TSGLEEGetType(TS ts,TSGLEEType *gleetype) 991 { 992 PetscErrorCode ierr; 993 994 PetscFunctionBegin; 995 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 996 ierr = PetscUseMethod(ts,"TSGLEEGetType_C",(TS,TSGLEEType*),(ts,gleetype));CHKERRQ(ierr); 997 PetscFunctionReturn(0); 998 } 999 1000 PetscErrorCode TSGLEEGetType_GLEE(TS ts,TSGLEEType *gleetype) 1001 { 1002 TS_GLEE *glee = (TS_GLEE*)ts->data; 1003 PetscErrorCode ierr; 1004 1005 PetscFunctionBegin; 1006 if (!glee->tableau) { 1007 ierr = TSGLEESetType(ts,TSGLEEDefaultType);CHKERRQ(ierr); 1008 } 1009 *gleetype = glee->tableau->name; 1010 PetscFunctionReturn(0); 1011 } 1012 PetscErrorCode TSGLEESetType_GLEE(TS ts,TSGLEEType gleetype) 1013 { 1014 TS_GLEE *glee = (TS_GLEE*)ts->data; 1015 PetscErrorCode ierr; 1016 PetscBool match; 1017 GLEETableauLink link; 1018 1019 PetscFunctionBegin; 1020 if (glee->tableau) { 1021 ierr = PetscStrcmp(glee->tableau->name,gleetype,&match);CHKERRQ(ierr); 1022 if (match) PetscFunctionReturn(0); 1023 } 1024 for (link = GLEETableauList; link; link=link->next) { 1025 ierr = PetscStrcmp(link->tab.name,gleetype,&match);CHKERRQ(ierr); 1026 if (match) { 1027 ierr = TSReset_GLEE(ts);CHKERRQ(ierr); 1028 glee->tableau = &link->tab; 1029 PetscFunctionReturn(0); 1030 } 1031 } 1032 SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_UNKNOWN_TYPE,"Could not find '%s'",gleetype); 1033 PetscFunctionReturn(0); 1034 } 1035 1036 static PetscErrorCode TSGetStages_GLEE(TS ts,PetscInt *ns,Vec **Y) 1037 { 1038 TS_GLEE *glee = (TS_GLEE*)ts->data; 1039 1040 PetscFunctionBegin; 1041 if (ns) *ns = glee->tableau->s; 1042 if (Y) *Y = glee->YStage; 1043 PetscFunctionReturn(0); 1044 } 1045 1046 PetscErrorCode TSGetSolutionComponents_GLEE(TS ts,PetscInt *n,Vec *Y) 1047 { 1048 TS_GLEE *glee = (TS_GLEE*)ts->data; 1049 GLEETableau tab = glee->tableau; 1050 PetscErrorCode ierr; 1051 1052 PetscFunctionBegin; 1053 if (!Y) *n = tab->r; 1054 else { 1055 if ((*n >= 0) && (*n < tab->r)) { 1056 ierr = VecCopy(glee->Y[*n],*Y); CHKERRQ(ierr); 1057 } else SETERRQ3(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_OUTOFRANGE,"Second argument (%d) out of range[%d,%d].",*n,0,tab->r-1); 1058 } 1059 PetscFunctionReturn(0); 1060 } 1061 1062 PetscErrorCode TSGetAuxSolution_GLEE(TS ts,Vec *X) 1063 { 1064 TS_GLEE *glee = (TS_GLEE*)ts->data; 1065 GLEETableau tab = glee->tableau; 1066 PetscReal *F = tab->Fembed; 1067 PetscInt r = tab->r; 1068 Vec *Y = glee->Y; 1069 PetscScalar *wr = glee->rwork; 1070 PetscInt i; 1071 PetscErrorCode ierr; 1072 1073 PetscFunctionBegin; 1074 ierr = VecZeroEntries(*X);CHKERRQ(ierr); 1075 for (i=0; i<r; i++) wr[i] = F[i]; 1076 ierr = VecMAXPY((*X),r,wr,Y);CHKERRQ(ierr); 1077 PetscFunctionReturn(0); 1078 } 1079 1080 PetscErrorCode TSGetTimeError_GLEE(TS ts,PetscInt n,Vec *X) 1081 { 1082 TS_GLEE *glee = (TS_GLEE*)ts->data; 1083 GLEETableau tab = glee->tableau; 1084 PetscReal *F = tab->Ferror; 1085 PetscInt r = tab->r; 1086 Vec *Y = glee->Y; 1087 PetscScalar *wr = glee->rwork; 1088 PetscInt i; 1089 PetscErrorCode ierr; 1090 1091 PetscFunctionBegin; 1092 ierr = VecZeroEntries(*X);CHKERRQ(ierr); 1093 if(n==0){ 1094 for (i=0; i<r; i++) wr[i] = F[i]; 1095 ierr = VecMAXPY((*X),r,wr,Y);CHKERRQ(ierr); 1096 } else if(n==-1) { 1097 *X=glee->yGErr; 1098 } 1099 PetscFunctionReturn(0); 1100 } 1101 1102 PetscErrorCode TSSetTimeError_GLEE(TS ts,Vec X) 1103 { 1104 TS_GLEE *glee = (TS_GLEE*)ts->data; 1105 GLEETableau tab = glee->tableau; 1106 PetscReal *S = tab->Serror; 1107 PetscInt r = tab->r,i; 1108 Vec *Y = glee->Y; 1109 PetscErrorCode ierr; 1110 1111 PetscFunctionBegin; 1112 if (r != 2) SETERRQ2(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"TSSetTimeError_GLEE not supported for '%s' with r=%D.",tab->name,tab->r); 1113 for (i=1; i<r; i++) { 1114 ierr = VecCopy(ts->vec_sol,Y[i]); CHKERRQ(ierr); 1115 ierr = VecAXPBY(Y[i],S[0],S[1],X); CHKERRQ(ierr); 1116 ierr = VecCopy(X,glee->yGErr); CHKERRQ(ierr); 1117 } 1118 PetscFunctionReturn(0); 1119 } 1120 1121 static PetscErrorCode TSDestroy_GLEE(TS ts) 1122 { 1123 PetscErrorCode ierr; 1124 1125 PetscFunctionBegin; 1126 ierr = TSReset_GLEE(ts);CHKERRQ(ierr); 1127 if (ts->dm) { 1128 ierr = DMCoarsenHookRemove(ts->dm,DMCoarsenHook_TSGLEE,DMRestrictHook_TSGLEE,ts);CHKERRQ(ierr); 1129 ierr = DMSubDomainHookRemove(ts->dm,DMSubDomainHook_TSGLEE,DMSubDomainRestrictHook_TSGLEE,ts);CHKERRQ(ierr); 1130 } 1131 ierr = PetscFree(ts->data);CHKERRQ(ierr); 1132 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSGLEEGetType_C",NULL);CHKERRQ(ierr); 1133 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSGLEESetType_C",NULL);CHKERRQ(ierr); 1134 PetscFunctionReturn(0); 1135 } 1136 1137 /* ------------------------------------------------------------ */ 1138 /*MC 1139 TSGLEE - ODE and DAE solver using General Linear with Error Estimation schemes 1140 1141 The user should provide the right hand side of the equation 1142 using TSSetRHSFunction(). 1143 1144 Notes: 1145 The default is TSGLEE35, it can be changed with TSGLEESetType() or -ts_glee_type 1146 1147 Level: beginner 1148 1149 .seealso: TSCreate(), TS, TSSetType(), TSGLEESetType(), TSGLEEGetType(), 1150 TSGLEE23, TTSGLEE24, TSGLEE35, TSGLEE25I, TSGLEEEXRK2A, 1151 TSGLEERK32G1, TSGLEERK285EX, TSGLEEType, TSGLEERegister() 1152 1153 M*/ 1154 PETSC_EXTERN PetscErrorCode TSCreate_GLEE(TS ts) 1155 { 1156 TS_GLEE *th; 1157 PetscErrorCode ierr; 1158 1159 PetscFunctionBegin; 1160 #if !defined(PETSC_USE_DYNAMIC_LIBRARIES) 1161 ierr = TSGLEEInitializePackage();CHKERRQ(ierr); 1162 #endif 1163 1164 ts->ops->reset = TSReset_GLEE; 1165 ts->ops->destroy = TSDestroy_GLEE; 1166 ts->ops->view = TSView_GLEE; 1167 ts->ops->load = TSLoad_GLEE; 1168 ts->ops->setup = TSSetUp_GLEE; 1169 ts->ops->step = TSStep_GLEE; 1170 ts->ops->interpolate = TSInterpolate_GLEE; 1171 ts->ops->evaluatestep = TSEvaluateStep_GLEE; 1172 ts->ops->setfromoptions = TSSetFromOptions_GLEE; 1173 ts->ops->getstages = TSGetStages_GLEE; 1174 ts->ops->snesfunction = SNESTSFormFunction_GLEE; 1175 ts->ops->snesjacobian = SNESTSFormJacobian_GLEE; 1176 ts->ops->getsolutioncomponents = TSGetSolutionComponents_GLEE; 1177 ts->ops->getauxsolution = TSGetAuxSolution_GLEE; 1178 ts->ops->gettimeerror = TSGetTimeError_GLEE; 1179 ts->ops->settimeerror = TSSetTimeError_GLEE; 1180 ts->ops->startingmethod = TSStartingMethod_GLEE; 1181 ts->default_adapt_type = TSADAPTGLEE; 1182 1183 ts->usessnes = PETSC_TRUE; 1184 1185 ierr = PetscNewLog(ts,&th);CHKERRQ(ierr); 1186 ts->data = (void*)th; 1187 1188 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSGLEEGetType_C",TSGLEEGetType_GLEE);CHKERRQ(ierr); 1189 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSGLEESetType_C",TSGLEESetType_GLEE);CHKERRQ(ierr); 1190 PetscFunctionReturn(0); 1191 } 1192