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