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 = glee->tableau; 854 PetscInt r=tab->r,i; 855 PetscReal *S=tab->S; 856 PetscErrorCode ierr; 857 858 PetscFunctionBegin; 859 for (i=0; i<r; i++) { 860 ierr = VecZeroEntries(glee->Y[i]);CHKERRQ(ierr); 861 ierr = VecAXPY(glee->Y[i],S[i],ts->vec_sol);CHKERRQ(ierr); 862 } 863 864 PetscFunctionReturn(0); 865 } 866 867 868 /*------------------------------------------------------------*/ 869 870 #undef __FUNCT__ 871 #define __FUNCT__ "TSSetFromOptions_GLEE" 872 static PetscErrorCode TSSetFromOptions_GLEE(PetscOptions *PetscOptionsObject,TS ts) 873 { 874 PetscErrorCode ierr; 875 char gleetype[256]; 876 877 PetscFunctionBegin; 878 ierr = PetscOptionsHead(PetscOptionsObject,"GLEE ODE solver options");CHKERRQ(ierr); 879 { 880 GLEETableauLink link; 881 PetscInt count,choice; 882 PetscBool flg; 883 const char **namelist; 884 885 ierr = PetscStrncpy(gleetype,TSGLEEDefaultType,sizeof(gleetype));CHKERRQ(ierr); 886 for (link=GLEETableauList,count=0; link; link=link->next,count++) ; 887 ierr = PetscMalloc1(count,&namelist);CHKERRQ(ierr); 888 for (link=GLEETableauList,count=0; link; link=link->next,count++) namelist[count] = link->tab.name; 889 ierr = PetscOptionsEList("-ts_glee_type","Family of GLEE method","TSGLEESetType",(const char*const*)namelist,count,gleetype,&choice,&flg);CHKERRQ(ierr); 890 ierr = TSGLEESetType(ts,flg ? namelist[choice] : gleetype);CHKERRQ(ierr); 891 ierr = PetscFree(namelist);CHKERRQ(ierr); 892 } 893 ierr = PetscOptionsTail();CHKERRQ(ierr); 894 PetscFunctionReturn(0); 895 } 896 897 #undef __FUNCT__ 898 #define __FUNCT__ "PetscFormatRealArray" 899 static PetscErrorCode PetscFormatRealArray(char buf[],size_t len,const char *fmt,PetscInt n,const PetscReal x[]) 900 { 901 PetscErrorCode ierr; 902 PetscInt i; 903 size_t left,count; 904 char *p; 905 906 PetscFunctionBegin; 907 for (i=0,p=buf,left=len; i<n; i++) { 908 ierr = PetscSNPrintfCount(p,left,fmt,&count,x[i]);CHKERRQ(ierr); 909 if (count >= left) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Insufficient space in buffer"); 910 left -= count; 911 p += count; 912 *p++ = ' '; 913 } 914 p[i ? 0 : -1] = 0; 915 PetscFunctionReturn(0); 916 } 917 918 #undef __FUNCT__ 919 #define __FUNCT__ "TSView_GLEE" 920 static PetscErrorCode TSView_GLEE(TS ts,PetscViewer viewer) 921 { 922 TS_GLEE *glee = (TS_GLEE*)ts->data; 923 GLEETableau tab = glee->tableau; 924 PetscBool iascii; 925 PetscErrorCode ierr; 926 TSAdapt adapt; 927 928 PetscFunctionBegin; 929 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 930 if (iascii) { 931 TSGLEEType gleetype; 932 char buf[512]; 933 ierr = TSGLEEGetType(ts,&gleetype);CHKERRQ(ierr); 934 ierr = PetscViewerASCIIPrintf(viewer," GLEE %s\n",gleetype);CHKERRQ(ierr); 935 ierr = PetscFormatRealArray(buf,sizeof(buf),"% 8.6f",tab->s,tab->c);CHKERRQ(ierr); 936 ierr = PetscViewerASCIIPrintf(viewer," Abscissa c = %s\n",buf);CHKERRQ(ierr); 937 /* Note: print out r as well */ 938 } 939 ierr = TSGetAdapt(ts,&adapt);CHKERRQ(ierr); 940 ierr = TSAdaptView(adapt,viewer);CHKERRQ(ierr); 941 PetscFunctionReturn(0); 942 } 943 944 #undef __FUNCT__ 945 #define __FUNCT__ "TSLoad_GLEE" 946 static PetscErrorCode TSLoad_GLEE(TS ts,PetscViewer viewer) 947 { 948 PetscErrorCode ierr; 949 SNES snes; 950 TSAdapt tsadapt; 951 952 PetscFunctionBegin; 953 ierr = TSGetAdapt(ts,&tsadapt);CHKERRQ(ierr); 954 ierr = TSAdaptLoad(tsadapt,viewer);CHKERRQ(ierr); 955 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 956 ierr = SNESLoad(snes,viewer);CHKERRQ(ierr); 957 /* function and Jacobian context for SNES when used with TS is always ts object */ 958 ierr = SNESSetFunction(snes,NULL,NULL,ts);CHKERRQ(ierr); 959 ierr = SNESSetJacobian(snes,NULL,NULL,NULL,ts);CHKERRQ(ierr); 960 PetscFunctionReturn(0); 961 } 962 963 #undef __FUNCT__ 964 #define __FUNCT__ "TSGLEESetType" 965 /*@C 966 TSGLEESetType - Set the type of GLEE scheme 967 968 Logically collective 969 970 Input Parameter: 971 + ts - timestepping context 972 - gleetype - type of GLEE-scheme 973 974 Level: intermediate 975 976 .seealso: TSGLEEGetType(), TSGLEE 977 @*/ 978 PetscErrorCode TSGLEESetType(TS ts,TSGLEEType gleetype) 979 { 980 PetscErrorCode ierr; 981 982 PetscFunctionBegin; 983 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 984 ierr = PetscTryMethod(ts,"TSGLEESetType_C",(TS,TSGLEEType),(ts,gleetype));CHKERRQ(ierr); 985 PetscFunctionReturn(0); 986 } 987 988 #undef __FUNCT__ 989 #define __FUNCT__ "TSGLEEGetType" 990 /*@C 991 TSGLEEGetType - Get the type of GLEE scheme 992 993 Logically collective 994 995 Input Parameter: 996 . ts - timestepping context 997 998 Output Parameter: 999 . gleetype - type of GLEE-scheme 1000 1001 Level: intermediate 1002 1003 .seealso: TSGLEESetType() 1004 @*/ 1005 PetscErrorCode TSGLEEGetType(TS ts,TSGLEEType *gleetype) 1006 { 1007 PetscErrorCode ierr; 1008 1009 PetscFunctionBegin; 1010 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1011 ierr = PetscUseMethod(ts,"TSGLEEGetType_C",(TS,TSGLEEType*),(ts,gleetype));CHKERRQ(ierr); 1012 PetscFunctionReturn(0); 1013 } 1014 1015 #undef __FUNCT__ 1016 #define __FUNCT__ "TSGLEEGetType_GLEE" 1017 PetscErrorCode TSGLEEGetType_GLEE(TS ts,TSGLEEType *gleetype) 1018 { 1019 TS_GLEE *glee = (TS_GLEE*)ts->data; 1020 PetscErrorCode ierr; 1021 1022 PetscFunctionBegin; 1023 if (!glee->tableau) { 1024 ierr = TSGLEESetType(ts,TSGLEEDefaultType);CHKERRQ(ierr); 1025 } 1026 *gleetype = glee->tableau->name; 1027 PetscFunctionReturn(0); 1028 } 1029 #undef __FUNCT__ 1030 #define __FUNCT__ "TSGLEESetType_GLEE" 1031 PetscErrorCode TSGLEESetType_GLEE(TS ts,TSGLEEType gleetype) 1032 { 1033 TS_GLEE *glee = (TS_GLEE*)ts->data; 1034 PetscErrorCode ierr; 1035 PetscBool match; 1036 GLEETableauLink link; 1037 1038 PetscFunctionBegin; 1039 if (glee->tableau) { 1040 ierr = PetscStrcmp(glee->tableau->name,gleetype,&match);CHKERRQ(ierr); 1041 if (match) PetscFunctionReturn(0); 1042 } 1043 for (link = GLEETableauList; link; link=link->next) { 1044 ierr = PetscStrcmp(link->tab.name,gleetype,&match);CHKERRQ(ierr); 1045 if (match) { 1046 ierr = TSReset_GLEE(ts);CHKERRQ(ierr); 1047 glee->tableau = &link->tab; 1048 PetscFunctionReturn(0); 1049 } 1050 } 1051 SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_UNKNOWN_TYPE,"Could not find '%s'",gleetype); 1052 PetscFunctionReturn(0); 1053 } 1054 1055 #undef __FUNCT__ 1056 #define __FUNCT__ "TSGetStages_GLEE" 1057 static PetscErrorCode TSGetStages_GLEE(TS ts,PetscInt *ns,Vec **Y) 1058 { 1059 TS_GLEE *glee = (TS_GLEE*)ts->data; 1060 1061 PetscFunctionBegin; 1062 *ns = glee->tableau->s; 1063 if(Y) *Y = glee->Y; 1064 PetscFunctionReturn(0); 1065 } 1066 1067 #undef __FUNCT__ 1068 #define __FUNCT__ "TSGetAuxSolution_GLEE" 1069 PetscErrorCode TSGetAuxSolution_GLEE(TS ts,PetscInt *n,Vec *Y) 1070 { 1071 TS_GLEE *glee = (TS_GLEE*)ts->data; 1072 GLEETableau tab = glee->tableau; 1073 PetscErrorCode ierr; 1074 1075 PetscFunctionBegin; 1076 if (!Y) *n = tab->r - 1; 1077 else { 1078 if ((*n >= 0) && (*n < tab->r-1)) { 1079 ierr = VecCopy(glee->Y[*n],*Y); CHKERRQ(ierr); 1080 } else { 1081 /* XXX: Error message for invalid value of n */ 1082 } 1083 } 1084 PetscFunctionReturn(0); 1085 } 1086 1087 /* ------------------------------------------------------------ */ 1088 /*MC 1089 TSGLEE - ODE and DAE solver using General Linear with Error Estimation schemes 1090 1091 The user should provide the right hand side of the equation 1092 using TSSetRHSFunction(). 1093 1094 Notes: 1095 The default is TSGLEE35, it can be changed with TSGLEESetType() or -ts_glee_type 1096 1097 Level: beginner 1098 1099 .seealso: TSCreate(), TS, TSSetType(), TSGLEESetType(), TSGLEEGetType(), 1100 TSGLEE23, TTSGLEE24, TSGLEE35, TSGLEE25I, TSGLEEEXRK2A, 1101 TSGLEERK32G1, TSGLEERK285EX, TSGLEEType, TSGLEERegister() 1102 1103 M*/ 1104 #undef __FUNCT__ 1105 #define __FUNCT__ "TSCreate_GLEE" 1106 PETSC_EXTERN PetscErrorCode TSCreate_GLEE(TS ts) 1107 { 1108 TS_GLEE *th; 1109 PetscErrorCode ierr; 1110 1111 PetscFunctionBegin; 1112 #if !defined(PETSC_USE_DYNAMIC_LIBRARIES) 1113 ierr = TSGLEEInitializePackage();CHKERRQ(ierr); 1114 #endif 1115 1116 ts->ops->reset = TSReset_GLEE; 1117 ts->ops->destroy = TSDestroy_GLEE; 1118 ts->ops->view = TSView_GLEE; 1119 ts->ops->load = TSLoad_GLEE; 1120 ts->ops->setup = TSSetUp_GLEE; 1121 ts->ops->step = TSStep_GLEE; 1122 ts->ops->interpolate = TSInterpolate_GLEE; 1123 ts->ops->evaluatestep = TSEvaluateStep_GLEE; 1124 ts->ops->setfromoptions = TSSetFromOptions_GLEE; 1125 ts->ops->getstages = TSGetStages_GLEE; 1126 ts->ops->snesfunction = SNESTSFormFunction_GLEE; 1127 ts->ops->snesjacobian = SNESTSFormJacobian_GLEE; 1128 ts->ops->getauxsolution = TSGetAuxSolution_GLEE; 1129 ts->ops->startingmethod = TSStartingMethod_GLEE; 1130 1131 ierr = PetscNewLog(ts,&th);CHKERRQ(ierr); 1132 ts->data = (void*)th; 1133 1134 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSGLEEGetType_C",TSGLEEGetType_GLEE);CHKERRQ(ierr); 1135 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSGLEESetType_C",TSGLEESetType_GLEE);CHKERRQ(ierr); 1136 PetscFunctionReturn(0); 1137 } 1138