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