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