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 PetscCall(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 PetscCall(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 PetscCall(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 PetscCall(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 PetscCall(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 PetscCall(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 PetscCall(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 PetscCall(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 GLEETableauLink link; 367 368 PetscFunctionBegin; 369 while ((link = GLEETableauList)) { 370 GLEETableau t = &link->tab; 371 GLEETableauList = link->next; 372 PetscCall(PetscFree5(t->A,t->B,t->U,t->V,t->c)); 373 PetscCall(PetscFree2(t->S,t->F)); 374 PetscCall(PetscFree (t->Fembed)); 375 PetscCall(PetscFree (t->Ferror)); 376 PetscCall(PetscFree (t->Serror)); 377 PetscCall(PetscFree (t->binterp)); 378 PetscCall(PetscFree (t->name)); 379 PetscCall(PetscFree (link)); 380 } 381 TSGLEERegisterAllCalled = PETSC_FALSE; 382 PetscFunctionReturn(0); 383 } 384 385 /*@C 386 TSGLEEInitializePackage - This function initializes everything in the TSGLEE package. It is called 387 from TSInitializePackage(). 388 389 Level: developer 390 391 .seealso: PetscInitialize() 392 @*/ 393 PetscErrorCode TSGLEEInitializePackage(void) 394 { 395 PetscFunctionBegin; 396 if (TSGLEEPackageInitialized) PetscFunctionReturn(0); 397 TSGLEEPackageInitialized = PETSC_TRUE; 398 PetscCall(TSGLEERegisterAll()); 399 PetscCall(PetscObjectComposedDataRegister(&explicit_stage_time_id)); 400 PetscCall(PetscRegisterFinalize(TSGLEEFinalizePackage)); 401 PetscFunctionReturn(0); 402 } 403 404 /*@C 405 TSGLEEFinalizePackage - This function destroys everything in the TSGLEE package. It is 406 called from PetscFinalize(). 407 408 Level: developer 409 410 .seealso: PetscFinalize() 411 @*/ 412 PetscErrorCode TSGLEEFinalizePackage(void) 413 { 414 PetscFunctionBegin; 415 TSGLEEPackageInitialized = PETSC_FALSE; 416 PetscCall(TSGLEERegisterDestroy()); 417 PetscFunctionReturn(0); 418 } 419 420 /*@C 421 TSGLEERegister - register an GLEE scheme by providing the entries in the Butcher tableau 422 423 Not Collective, but the same schemes should be registered on all processes on which they will be used 424 425 Input Parameters: 426 + name - identifier for method 427 . order - order of method 428 . s - number of stages 429 . r - number of steps 430 . gamma - LTE ratio 431 . A - stage coefficients (dimension s*s, row-major) 432 . B - step completion coefficients (dimension r*s, row-major) 433 . U - method coefficients (dimension s*r, row-major) 434 . V - method coefficients (dimension r*r, row-major) 435 . S - starting coefficients 436 . F - finishing coefficients 437 . c - abscissa (dimension s; NULL to use row sums of A) 438 . Fembed - step completion coefficients for embedded method 439 . Ferror - error computation coefficients 440 . Serror - error initialization coefficients 441 . pinterp - order of interpolation (0 if unavailable) 442 - binterp - array of interpolation coefficients (NULL if unavailable) 443 444 Notes: 445 Several GLEE methods are provided, this function is only needed to create new methods. 446 447 Level: advanced 448 449 .seealso: TSGLEE 450 @*/ 451 PetscErrorCode TSGLEERegister(TSGLEEType name,PetscInt order,PetscInt s, PetscInt r, 452 PetscReal gamma, 453 const PetscReal A[],const PetscReal B[], 454 const PetscReal U[],const PetscReal V[], 455 const PetscReal S[],const PetscReal F[], 456 const PetscReal c[], 457 const PetscReal Fembed[],const PetscReal Ferror[], 458 const PetscReal Serror[], 459 PetscInt pinterp, const PetscReal binterp[]) 460 { 461 GLEETableauLink link; 462 GLEETableau t; 463 PetscInt i,j; 464 465 PetscFunctionBegin; 466 PetscCall(TSGLEEInitializePackage()); 467 PetscCall(PetscNew(&link)); 468 t = &link->tab; 469 PetscCall(PetscStrallocpy(name,&t->name)); 470 t->order = order; 471 t->s = s; 472 t->r = r; 473 t->gamma = gamma; 474 PetscCall(PetscMalloc5(s*s,&t->A,r*r,&t->V,s,&t->c,r*s,&t->B,s*r,&t->U)); 475 PetscCall(PetscMalloc2(r,&t->S,r,&t->F)); 476 PetscCall(PetscArraycpy(t->A,A,s*s)); 477 PetscCall(PetscArraycpy(t->B,B,r*s)); 478 PetscCall(PetscArraycpy(t->U,U,s*r)); 479 PetscCall(PetscArraycpy(t->V,V,r*r)); 480 PetscCall(PetscArraycpy(t->S,S,r)); 481 PetscCall(PetscArraycpy(t->F,F,r)); 482 if (c) { 483 PetscCall(PetscArraycpy(t->c,c,s)); 484 } else { 485 for (i=0; i<s; i++) for (j=0,t->c[i]=0; j<s; j++) t->c[i] += A[i*s+j]; 486 } 487 PetscCall(PetscMalloc1(r,&t->Fembed)); 488 PetscCall(PetscMalloc1(r,&t->Ferror)); 489 PetscCall(PetscMalloc1(r,&t->Serror)); 490 PetscCall(PetscArraycpy(t->Fembed,Fembed,r)); 491 PetscCall(PetscArraycpy(t->Ferror,Ferror,r)); 492 PetscCall(PetscArraycpy(t->Serror,Serror,r)); 493 t->pinterp = pinterp; 494 PetscCall(PetscMalloc1(s*pinterp,&t->binterp)); 495 PetscCall(PetscArraycpy(t->binterp,binterp,s*pinterp)); 496 497 link->next = GLEETableauList; 498 GLEETableauList = link; 499 PetscFunctionReturn(0); 500 } 501 502 static PetscErrorCode TSEvaluateStep_GLEE(TS ts,PetscInt order,Vec X,PetscBool *done) 503 { 504 TS_GLEE *glee = (TS_GLEE*) ts->data; 505 GLEETableau tab = glee->tableau; 506 PetscReal h, *B = tab->B, *V = tab->V, 507 *F = tab->F, 508 *Fembed = tab->Fembed; 509 PetscInt s = tab->s, r = tab->r, i, j; 510 Vec *Y = glee->Y, *YdotStage = glee->YdotStage; 511 PetscScalar *ws = glee->swork, *wr = glee->rwork; 512 513 PetscFunctionBegin; 514 515 switch (glee->status) { 516 case TS_STEP_INCOMPLETE: 517 case TS_STEP_PENDING: 518 h = ts->time_step; break; 519 case TS_STEP_COMPLETE: 520 h = ts->ptime - ts->ptime_prev; break; 521 default: SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_PLIB,"Invalid TSStepStatus"); 522 } 523 524 if (order == tab->order) { 525 526 /* Note: Irrespective of whether status is TS_STEP_INCOMPLETE 527 or TS_STEP_COMPLETE, glee->X has the solution at the 528 beginning of the time step. So no need to roll-back. 529 */ 530 if (glee->status == TS_STEP_INCOMPLETE) { 531 for (i=0; i<r; i++) { 532 PetscCall(VecZeroEntries(Y[i])); 533 for (j=0; j<r; j++) wr[j] = V[i*r+j]; 534 PetscCall(VecMAXPY(Y[i],r,wr,glee->X)); 535 for (j=0; j<s; j++) ws[j] = h*B[i*s+j]; 536 PetscCall(VecMAXPY(Y[i],s,ws,YdotStage)); 537 } 538 PetscCall(VecZeroEntries(X)); 539 for (j=0; j<r; j++) wr[j] = F[j]; 540 PetscCall(VecMAXPY(X,r,wr,Y)); 541 } else PetscCall(VecCopy(ts->vec_sol,X)); 542 PetscFunctionReturn(0); 543 544 } else if (order == tab->order-1) { 545 546 /* Complete with the embedded method (Fembed) */ 547 for (i=0; i<r; i++) { 548 PetscCall(VecZeroEntries(Y[i])); 549 for (j=0; j<r; j++) wr[j] = V[i*r+j]; 550 PetscCall(VecMAXPY(Y[i],r,wr,glee->X)); 551 for (j=0; j<s; j++) ws[j] = h*B[i*s+j]; 552 PetscCall(VecMAXPY(Y[i],s,ws,YdotStage)); 553 } 554 PetscCall(VecZeroEntries(X)); 555 for (j=0; j<r; j++) wr[j] = Fembed[j]; 556 PetscCall(VecMAXPY(X,r,wr,Y)); 557 558 if (done) *done = PETSC_TRUE; 559 PetscFunctionReturn(0); 560 } 561 if (done) *done = PETSC_FALSE; 562 else SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"GLEE '%s' of order %" PetscInt_FMT " cannot evaluate step at order %" PetscInt_FMT,tab->name,tab->order,order); 563 PetscFunctionReturn(0); 564 } 565 566 static PetscErrorCode TSStep_GLEE(TS ts) 567 { 568 TS_GLEE *glee = (TS_GLEE*)ts->data; 569 GLEETableau tab = glee->tableau; 570 const PetscInt s = tab->s, r = tab->r; 571 PetscReal *A = tab->A, *U = tab->U, 572 *F = tab->F, 573 *c = tab->c; 574 Vec *Y = glee->Y, *X = glee->X, 575 *YStage = glee->YStage, 576 *YdotStage = glee->YdotStage, 577 W = glee->W; 578 SNES snes; 579 PetscScalar *ws = glee->swork, *wr = glee->rwork; 580 TSAdapt adapt; 581 PetscInt i,j,reject,next_scheme,its,lits; 582 PetscReal next_time_step; 583 PetscReal t; 584 PetscBool accept; 585 586 PetscFunctionBegin; 587 PetscCall(PetscCitationsRegister(citation,&cited)); 588 589 for (i=0; i<r; i++) PetscCall(VecCopy(Y[i],X[i])); 590 591 PetscCall(TSGetSNES(ts,&snes)); 592 next_time_step = ts->time_step; 593 t = ts->ptime; 594 accept = PETSC_TRUE; 595 glee->status = TS_STEP_INCOMPLETE; 596 597 for (reject=0; reject<ts->max_reject && !ts->reason; reject++,ts->reject++) { 598 599 PetscReal h = ts->time_step; 600 PetscCall(TSPreStep(ts)); 601 602 for (i=0; i<s; i++) { 603 604 glee->stage_time = t + h*c[i]; 605 PetscCall(TSPreStage(ts,glee->stage_time)); 606 607 if (A[i*s+i] == 0) { /* Explicit stage */ 608 PetscCall(VecZeroEntries(YStage[i])); 609 for (j=0; j<r; j++) wr[j] = U[i*r+j]; 610 PetscCall(VecMAXPY(YStage[i],r,wr,X)); 611 for (j=0; j<i; j++) ws[j] = h*A[i*s+j]; 612 PetscCall(VecMAXPY(YStage[i],i,ws,YdotStage)); 613 } else { /* Implicit stage */ 614 glee->scoeff = 1.0/A[i*s+i]; 615 /* compute right-hand-side */ 616 PetscCall(VecZeroEntries(W)); 617 for (j=0; j<r; j++) wr[j] = U[i*r+j]; 618 PetscCall(VecMAXPY(W,r,wr,X)); 619 for (j=0; j<i; j++) ws[j] = h*A[i*s+j]; 620 PetscCall(VecMAXPY(W,i,ws,YdotStage)); 621 PetscCall(VecScale(W,glee->scoeff/h)); 622 /* set initial guess */ 623 PetscCall(VecCopy(i>0 ? YStage[i-1] : ts->vec_sol,YStage[i])); 624 /* solve for this stage */ 625 PetscCall(SNESSolve(snes,W,YStage[i])); 626 PetscCall(SNESGetIterationNumber(snes,&its)); 627 PetscCall(SNESGetLinearSolveIterations(snes,&lits)); 628 ts->snes_its += its; ts->ksp_its += lits; 629 } 630 PetscCall(TSGetAdapt(ts,&adapt)); 631 PetscCall(TSAdaptCheckStage(adapt,ts,glee->stage_time,YStage[i],&accept)); 632 if (!accept) goto reject_step; 633 PetscCall(TSPostStage(ts,glee->stage_time,i,YStage)); 634 PetscCall(TSComputeRHSFunction(ts,t+h*c[i],YStage[i],YdotStage[i])); 635 } 636 PetscCall(TSEvaluateStep(ts,tab->order,ts->vec_sol,NULL)); 637 glee->status = TS_STEP_PENDING; 638 639 /* Register only the current method as a candidate because we're not supporting multiple candidates yet. */ 640 PetscCall(TSGetAdapt(ts,&adapt)); 641 PetscCall(TSAdaptCandidatesClear(adapt)); 642 PetscCall(TSAdaptCandidateAdd(adapt,tab->name,tab->order,1,tab->ccfl,(PetscReal)tab->s,PETSC_TRUE)); 643 PetscCall(TSAdaptChoose(adapt,ts,ts->time_step,&next_scheme,&next_time_step,&accept)); 644 if (accept) { 645 /* ignore next_scheme for now */ 646 ts->ptime += ts->time_step; 647 ts->time_step = next_time_step; 648 glee->status = TS_STEP_COMPLETE; 649 /* compute and store the global error */ 650 /* Note: this is not needed if TSAdaptGLEE is not used */ 651 PetscCall(TSGetTimeError(ts,0,&(glee->yGErr))); 652 PetscCall(PetscObjectComposedDataSetReal((PetscObject)ts->vec_sol,explicit_stage_time_id,ts->ptime)); 653 break; 654 } else { /* Roll back the current step */ 655 for (j=0; j<r; j++) wr[j] = F[j]; 656 PetscCall(VecMAXPY(ts->vec_sol,r,wr,X)); 657 ts->time_step = next_time_step; 658 glee->status = TS_STEP_INCOMPLETE; 659 } 660 reject_step: continue; 661 } 662 if (glee->status != TS_STEP_COMPLETE && !ts->reason) ts->reason = TS_DIVERGED_STEP_REJECTED; 663 PetscFunctionReturn(0); 664 } 665 666 static PetscErrorCode TSInterpolate_GLEE(TS ts,PetscReal itime,Vec X) 667 { 668 TS_GLEE *glee = (TS_GLEE*)ts->data; 669 PetscInt s=glee->tableau->s, pinterp=glee->tableau->pinterp,i,j; 670 PetscReal h,tt,t; 671 PetscScalar *b; 672 const PetscReal *B = glee->tableau->binterp; 673 674 PetscFunctionBegin; 675 PetscCheck(B,PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"TSGLEE %s does not have an interpolation formula",glee->tableau->name); 676 switch (glee->status) { 677 case TS_STEP_INCOMPLETE: 678 case TS_STEP_PENDING: 679 h = ts->time_step; 680 t = (itime - ts->ptime)/h; 681 break; 682 case TS_STEP_COMPLETE: 683 h = ts->ptime - ts->ptime_prev; 684 t = (itime - ts->ptime)/h + 1; /* In the interval [0,1] */ 685 break; 686 default: SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_PLIB,"Invalid TSStepStatus"); 687 } 688 PetscCall(PetscMalloc1(s,&b)); 689 for (i=0; i<s; i++) b[i] = 0; 690 for (j=0,tt=t; j<pinterp; j++,tt*=t) { 691 for (i=0; i<s; i++) { 692 b[i] += h * B[i*pinterp+j] * tt; 693 } 694 } 695 PetscCall(VecCopy(glee->YStage[0],X)); 696 PetscCall(VecMAXPY(X,s,b,glee->YdotStage)); 697 PetscCall(PetscFree(b)); 698 PetscFunctionReturn(0); 699 } 700 701 /*------------------------------------------------------------*/ 702 static PetscErrorCode TSReset_GLEE(TS ts) 703 { 704 TS_GLEE *glee = (TS_GLEE*)ts->data; 705 PetscInt s, r; 706 707 PetscFunctionBegin; 708 if (!glee->tableau) PetscFunctionReturn(0); 709 s = glee->tableau->s; 710 r = glee->tableau->r; 711 PetscCall(VecDestroyVecs(r,&glee->Y)); 712 PetscCall(VecDestroyVecs(r,&glee->X)); 713 PetscCall(VecDestroyVecs(s,&glee->YStage)); 714 PetscCall(VecDestroyVecs(s,&glee->YdotStage)); 715 PetscCall(VecDestroy(&glee->Ydot)); 716 PetscCall(VecDestroy(&glee->yGErr)); 717 PetscCall(VecDestroy(&glee->W)); 718 PetscCall(PetscFree2(glee->swork,glee->rwork)); 719 PetscFunctionReturn(0); 720 } 721 722 static PetscErrorCode TSGLEEGetVecs(TS ts,DM dm,Vec *Ydot) 723 { 724 TS_GLEE *glee = (TS_GLEE*)ts->data; 725 726 PetscFunctionBegin; 727 if (Ydot) { 728 if (dm && dm != ts->dm) { 729 PetscCall(DMGetNamedGlobalVector(dm,"TSGLEE_Ydot",Ydot)); 730 } else *Ydot = glee->Ydot; 731 } 732 PetscFunctionReturn(0); 733 } 734 735 static PetscErrorCode TSGLEERestoreVecs(TS ts,DM dm,Vec *Ydot) 736 { 737 PetscFunctionBegin; 738 if (Ydot) { 739 if (dm && dm != ts->dm) { 740 PetscCall(DMRestoreNamedGlobalVector(dm,"TSGLEE_Ydot",Ydot)); 741 } 742 } 743 PetscFunctionReturn(0); 744 } 745 746 /* 747 This defines the nonlinear equation that is to be solved with SNES 748 */ 749 static PetscErrorCode SNESTSFormFunction_GLEE(SNES snes,Vec X,Vec F,TS ts) 750 { 751 TS_GLEE *glee = (TS_GLEE*)ts->data; 752 DM dm,dmsave; 753 Vec Ydot; 754 PetscReal shift = glee->scoeff / ts->time_step; 755 756 PetscFunctionBegin; 757 PetscCall(SNESGetDM(snes,&dm)); 758 PetscCall(TSGLEEGetVecs(ts,dm,&Ydot)); 759 /* Set Ydot = shift*X */ 760 PetscCall(VecCopy(X,Ydot)); 761 PetscCall(VecScale(Ydot,shift)); 762 dmsave = ts->dm; 763 ts->dm = dm; 764 765 PetscCall(TSComputeIFunction(ts,glee->stage_time,X,Ydot,F,PETSC_FALSE)); 766 767 ts->dm = dmsave; 768 PetscCall(TSGLEERestoreVecs(ts,dm,&Ydot)); 769 PetscFunctionReturn(0); 770 } 771 772 static PetscErrorCode SNESTSFormJacobian_GLEE(SNES snes,Vec X,Mat A,Mat B,TS ts) 773 { 774 TS_GLEE *glee = (TS_GLEE*)ts->data; 775 DM dm,dmsave; 776 Vec Ydot; 777 PetscReal shift = glee->scoeff / ts->time_step; 778 779 PetscFunctionBegin; 780 PetscCall(SNESGetDM(snes,&dm)); 781 PetscCall(TSGLEEGetVecs(ts,dm,&Ydot)); 782 /* glee->Ydot has already been computed in SNESTSFormFunction_GLEE (SNES guarantees this) */ 783 dmsave = ts->dm; 784 ts->dm = dm; 785 786 PetscCall(TSComputeIJacobian(ts,glee->stage_time,X,Ydot,shift,A,B,PETSC_FALSE)); 787 788 ts->dm = dmsave; 789 PetscCall(TSGLEERestoreVecs(ts,dm,&Ydot)); 790 PetscFunctionReturn(0); 791 } 792 793 static PetscErrorCode DMCoarsenHook_TSGLEE(DM fine,DM coarse,void *ctx) 794 { 795 PetscFunctionBegin; 796 PetscFunctionReturn(0); 797 } 798 799 static PetscErrorCode DMRestrictHook_TSGLEE(DM fine,Mat restrct,Vec rscale,Mat inject,DM coarse,void *ctx) 800 { 801 PetscFunctionBegin; 802 PetscFunctionReturn(0); 803 } 804 805 static PetscErrorCode DMSubDomainHook_TSGLEE(DM dm,DM subdm,void *ctx) 806 { 807 PetscFunctionBegin; 808 PetscFunctionReturn(0); 809 } 810 811 static PetscErrorCode DMSubDomainRestrictHook_TSGLEE(DM dm,VecScatter gscat,VecScatter lscat,DM subdm,void *ctx) 812 { 813 PetscFunctionBegin; 814 PetscFunctionReturn(0); 815 } 816 817 static PetscErrorCode TSSetUp_GLEE(TS ts) 818 { 819 TS_GLEE *glee = (TS_GLEE*)ts->data; 820 GLEETableau tab; 821 PetscInt s,r; 822 DM dm; 823 824 PetscFunctionBegin; 825 if (!glee->tableau) { 826 PetscCall(TSGLEESetType(ts,TSGLEEDefaultType)); 827 } 828 tab = glee->tableau; 829 s = tab->s; 830 r = tab->r; 831 PetscCall(VecDuplicateVecs(ts->vec_sol,r,&glee->Y)); 832 PetscCall(VecDuplicateVecs(ts->vec_sol,r,&glee->X)); 833 PetscCall(VecDuplicateVecs(ts->vec_sol,s,&glee->YStage)); 834 PetscCall(VecDuplicateVecs(ts->vec_sol,s,&glee->YdotStage)); 835 PetscCall(VecDuplicate(ts->vec_sol,&glee->Ydot)); 836 PetscCall(VecDuplicate(ts->vec_sol,&glee->yGErr)); 837 PetscCall(VecZeroEntries(glee->yGErr)); 838 PetscCall(VecDuplicate(ts->vec_sol,&glee->W)); 839 PetscCall(PetscMalloc2(s,&glee->swork,r,&glee->rwork)); 840 PetscCall(TSGetDM(ts,&dm)); 841 PetscCall(DMCoarsenHookAdd(dm,DMCoarsenHook_TSGLEE,DMRestrictHook_TSGLEE,ts)); 842 PetscCall(DMSubDomainHookAdd(dm,DMSubDomainHook_TSGLEE,DMSubDomainRestrictHook_TSGLEE,ts)); 843 PetscFunctionReturn(0); 844 } 845 846 PetscErrorCode TSStartingMethod_GLEE(TS ts) 847 { 848 TS_GLEE *glee = (TS_GLEE*)ts->data; 849 GLEETableau tab = glee->tableau; 850 PetscInt r=tab->r,i; 851 PetscReal *S=tab->S; 852 853 PetscFunctionBegin; 854 for (i=0; i<r; i++) { 855 PetscCall(VecZeroEntries(glee->Y[i])); 856 PetscCall(VecAXPY(glee->Y[i],S[i],ts->vec_sol)); 857 } 858 859 PetscFunctionReturn(0); 860 } 861 862 /*------------------------------------------------------------*/ 863 864 static PetscErrorCode TSSetFromOptions_GLEE(PetscOptionItems *PetscOptionsObject,TS ts) 865 { 866 char gleetype[256]; 867 868 PetscFunctionBegin; 869 PetscOptionsHeadBegin(PetscOptionsObject,"GLEE ODE solver options"); 870 { 871 GLEETableauLink link; 872 PetscInt count,choice; 873 PetscBool flg; 874 const char **namelist; 875 876 PetscCall(PetscStrncpy(gleetype,TSGLEEDefaultType,sizeof(gleetype))); 877 for (link=GLEETableauList,count=0; link; link=link->next,count++) ; 878 PetscCall(PetscMalloc1(count,(char***)&namelist)); 879 for (link=GLEETableauList,count=0; link; link=link->next,count++) namelist[count] = link->tab.name; 880 PetscCall(PetscOptionsEList("-ts_glee_type","Family of GLEE method","TSGLEESetType",(const char*const*)namelist,count,gleetype,&choice,&flg)); 881 PetscCall(TSGLEESetType(ts,flg ? namelist[choice] : gleetype)); 882 PetscCall(PetscFree(namelist)); 883 } 884 PetscOptionsHeadEnd(); 885 PetscFunctionReturn(0); 886 } 887 888 static PetscErrorCode TSView_GLEE(TS ts,PetscViewer viewer) 889 { 890 TS_GLEE *glee = (TS_GLEE*)ts->data; 891 GLEETableau tab = glee->tableau; 892 PetscBool iascii; 893 894 PetscFunctionBegin; 895 PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii)); 896 if (iascii) { 897 TSGLEEType gleetype; 898 char buf[512]; 899 PetscCall(TSGLEEGetType(ts,&gleetype)); 900 PetscCall(PetscViewerASCIIPrintf(viewer," GLEE type %s\n",gleetype)); 901 PetscCall(PetscFormatRealArray(buf,sizeof(buf),"% 8.6f",tab->s,tab->c)); 902 PetscCall(PetscViewerASCIIPrintf(viewer," Abscissa c = %s\n",buf)); 903 /* Note: print out r as well */ 904 } 905 PetscFunctionReturn(0); 906 } 907 908 static PetscErrorCode TSLoad_GLEE(TS ts,PetscViewer viewer) 909 { 910 SNES snes; 911 TSAdapt tsadapt; 912 913 PetscFunctionBegin; 914 PetscCall(TSGetAdapt(ts,&tsadapt)); 915 PetscCall(TSAdaptLoad(tsadapt,viewer)); 916 PetscCall(TSGetSNES(ts,&snes)); 917 PetscCall(SNESLoad(snes,viewer)); 918 /* function and Jacobian context for SNES when used with TS is always ts object */ 919 PetscCall(SNESSetFunction(snes,NULL,NULL,ts)); 920 PetscCall(SNESSetJacobian(snes,NULL,NULL,NULL,ts)); 921 PetscFunctionReturn(0); 922 } 923 924 /*@C 925 TSGLEESetType - Set the type of GLEE scheme 926 927 Logically collective 928 929 Input Parameters: 930 + ts - timestepping context 931 - gleetype - type of GLEE-scheme 932 933 Level: intermediate 934 935 .seealso: TSGLEEGetType(), TSGLEE 936 @*/ 937 PetscErrorCode TSGLEESetType(TS ts,TSGLEEType gleetype) 938 { 939 PetscFunctionBegin; 940 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 941 PetscValidCharPointer(gleetype,2); 942 PetscTryMethod(ts,"TSGLEESetType_C",(TS,TSGLEEType),(ts,gleetype)); 943 PetscFunctionReturn(0); 944 } 945 946 /*@C 947 TSGLEEGetType - Get the type of GLEE scheme 948 949 Logically collective 950 951 Input Parameter: 952 . ts - timestepping context 953 954 Output Parameter: 955 . gleetype - type of GLEE-scheme 956 957 Level: intermediate 958 959 .seealso: TSGLEESetType() 960 @*/ 961 PetscErrorCode TSGLEEGetType(TS ts,TSGLEEType *gleetype) 962 { 963 PetscFunctionBegin; 964 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 965 PetscUseMethod(ts,"TSGLEEGetType_C",(TS,TSGLEEType*),(ts,gleetype)); 966 PetscFunctionReturn(0); 967 } 968 969 PetscErrorCode TSGLEEGetType_GLEE(TS ts,TSGLEEType *gleetype) 970 { 971 TS_GLEE *glee = (TS_GLEE*)ts->data; 972 973 PetscFunctionBegin; 974 if (!glee->tableau) { 975 PetscCall(TSGLEESetType(ts,TSGLEEDefaultType)); 976 } 977 *gleetype = glee->tableau->name; 978 PetscFunctionReturn(0); 979 } 980 PetscErrorCode TSGLEESetType_GLEE(TS ts,TSGLEEType gleetype) 981 { 982 TS_GLEE *glee = (TS_GLEE*)ts->data; 983 PetscBool match; 984 GLEETableauLink link; 985 986 PetscFunctionBegin; 987 if (glee->tableau) { 988 PetscCall(PetscStrcmp(glee->tableau->name,gleetype,&match)); 989 if (match) PetscFunctionReturn(0); 990 } 991 for (link = GLEETableauList; link; link=link->next) { 992 PetscCall(PetscStrcmp(link->tab.name,gleetype,&match)); 993 if (match) { 994 PetscCall(TSReset_GLEE(ts)); 995 glee->tableau = &link->tab; 996 PetscFunctionReturn(0); 997 } 998 } 999 SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_UNKNOWN_TYPE,"Could not find '%s'",gleetype); 1000 } 1001 1002 static PetscErrorCode TSGetStages_GLEE(TS ts,PetscInt *ns,Vec **Y) 1003 { 1004 TS_GLEE *glee = (TS_GLEE*)ts->data; 1005 1006 PetscFunctionBegin; 1007 if (ns) *ns = glee->tableau->s; 1008 if (Y) *Y = glee->YStage; 1009 PetscFunctionReturn(0); 1010 } 1011 1012 PetscErrorCode TSGetSolutionComponents_GLEE(TS ts,PetscInt *n,Vec *Y) 1013 { 1014 TS_GLEE *glee = (TS_GLEE*)ts->data; 1015 GLEETableau tab = glee->tableau; 1016 1017 PetscFunctionBegin; 1018 if (!Y) *n = tab->r; 1019 else { 1020 if ((*n >= 0) && (*n < tab->r)) { 1021 PetscCall(VecCopy(glee->Y[*n],*Y)); 1022 } else SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_OUTOFRANGE,"Second argument (%" PetscInt_FMT ") out of range[0,%" PetscInt_FMT "].",*n,tab->r-1); 1023 } 1024 PetscFunctionReturn(0); 1025 } 1026 1027 PetscErrorCode TSGetAuxSolution_GLEE(TS ts,Vec *X) 1028 { 1029 TS_GLEE *glee = (TS_GLEE*)ts->data; 1030 GLEETableau tab = glee->tableau; 1031 PetscReal *F = tab->Fembed; 1032 PetscInt r = tab->r; 1033 Vec *Y = glee->Y; 1034 PetscScalar *wr = glee->rwork; 1035 PetscInt i; 1036 1037 PetscFunctionBegin; 1038 PetscCall(VecZeroEntries(*X)); 1039 for (i=0; i<r; i++) wr[i] = F[i]; 1040 PetscCall(VecMAXPY((*X),r,wr,Y)); 1041 PetscFunctionReturn(0); 1042 } 1043 1044 PetscErrorCode TSGetTimeError_GLEE(TS ts,PetscInt n,Vec *X) 1045 { 1046 TS_GLEE *glee = (TS_GLEE*)ts->data; 1047 GLEETableau tab = glee->tableau; 1048 PetscReal *F = tab->Ferror; 1049 PetscInt r = tab->r; 1050 Vec *Y = glee->Y; 1051 PetscScalar *wr = glee->rwork; 1052 PetscInt i; 1053 1054 PetscFunctionBegin; 1055 PetscCall(VecZeroEntries(*X)); 1056 if (n==0) { 1057 for (i=0; i<r; i++) wr[i] = F[i]; 1058 PetscCall(VecMAXPY((*X),r,wr,Y)); 1059 } else if (n==-1) { 1060 *X=glee->yGErr; 1061 } 1062 PetscFunctionReturn(0); 1063 } 1064 1065 PetscErrorCode TSSetTimeError_GLEE(TS ts,Vec X) 1066 { 1067 TS_GLEE *glee = (TS_GLEE*)ts->data; 1068 GLEETableau tab = glee->tableau; 1069 PetscReal *S = tab->Serror; 1070 PetscInt r = tab->r,i; 1071 Vec *Y = glee->Y; 1072 1073 PetscFunctionBegin; 1074 PetscCheck(r == 2,PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"TSSetTimeError_GLEE not supported for '%s' with r=%" PetscInt_FMT ".",tab->name,tab->r); 1075 for (i=1; i<r; i++) { 1076 PetscCall(VecCopy(ts->vec_sol,Y[i])); 1077 PetscCall(VecAXPBY(Y[i],S[0],S[1],X)); 1078 PetscCall(VecCopy(X,glee->yGErr)); 1079 } 1080 PetscFunctionReturn(0); 1081 } 1082 1083 static PetscErrorCode TSDestroy_GLEE(TS ts) 1084 { 1085 PetscFunctionBegin; 1086 PetscCall(TSReset_GLEE(ts)); 1087 if (ts->dm) { 1088 PetscCall(DMCoarsenHookRemove(ts->dm,DMCoarsenHook_TSGLEE,DMRestrictHook_TSGLEE,ts)); 1089 PetscCall(DMSubDomainHookRemove(ts->dm,DMSubDomainHook_TSGLEE,DMSubDomainRestrictHook_TSGLEE,ts)); 1090 } 1091 PetscCall(PetscFree(ts->data)); 1092 PetscCall(PetscObjectComposeFunction((PetscObject)ts,"TSGLEEGetType_C",NULL)); 1093 PetscCall(PetscObjectComposeFunction((PetscObject)ts,"TSGLEESetType_C",NULL)); 1094 PetscFunctionReturn(0); 1095 } 1096 1097 /* ------------------------------------------------------------ */ 1098 /*MC 1099 TSGLEE - ODE and DAE solver using General Linear with Error Estimation schemes 1100 1101 The user should provide the right hand side of the equation 1102 using TSSetRHSFunction(). 1103 1104 Notes: 1105 The default is TSGLEE35, it can be changed with TSGLEESetType() or -ts_glee_type 1106 1107 Level: beginner 1108 1109 .seealso: TSCreate(), TS, TSSetType(), TSGLEESetType(), TSGLEEGetType(), 1110 TSGLEE23, TTSGLEE24, TSGLEE35, TSGLEE25I, TSGLEEEXRK2A, 1111 TSGLEERK32G1, TSGLEERK285EX, TSGLEEType, TSGLEERegister() 1112 1113 M*/ 1114 PETSC_EXTERN PetscErrorCode TSCreate_GLEE(TS ts) 1115 { 1116 TS_GLEE *th; 1117 1118 PetscFunctionBegin; 1119 PetscCall(TSGLEEInitializePackage()); 1120 1121 ts->ops->reset = TSReset_GLEE; 1122 ts->ops->destroy = TSDestroy_GLEE; 1123 ts->ops->view = TSView_GLEE; 1124 ts->ops->load = TSLoad_GLEE; 1125 ts->ops->setup = TSSetUp_GLEE; 1126 ts->ops->step = TSStep_GLEE; 1127 ts->ops->interpolate = TSInterpolate_GLEE; 1128 ts->ops->evaluatestep = TSEvaluateStep_GLEE; 1129 ts->ops->setfromoptions = TSSetFromOptions_GLEE; 1130 ts->ops->getstages = TSGetStages_GLEE; 1131 ts->ops->snesfunction = SNESTSFormFunction_GLEE; 1132 ts->ops->snesjacobian = SNESTSFormJacobian_GLEE; 1133 ts->ops->getsolutioncomponents = TSGetSolutionComponents_GLEE; 1134 ts->ops->getauxsolution = TSGetAuxSolution_GLEE; 1135 ts->ops->gettimeerror = TSGetTimeError_GLEE; 1136 ts->ops->settimeerror = TSSetTimeError_GLEE; 1137 ts->ops->startingmethod = TSStartingMethod_GLEE; 1138 ts->default_adapt_type = TSADAPTGLEE; 1139 1140 ts->usessnes = PETSC_TRUE; 1141 1142 PetscCall(PetscNewLog(ts,&th)); 1143 ts->data = (void*)th; 1144 1145 PetscCall(PetscObjectComposeFunction((PetscObject)ts,"TSGLEEGetType_C",TSGLEEGetType_GLEE)); 1146 PetscCall(PetscObjectComposeFunction((PetscObject)ts,"TSGLEESetType_C",TSGLEESetType_GLEE)); 1147 PetscFunctionReturn(0); 1148 } 1149