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