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