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