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