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