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 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 if (done) *done = PETSC_FALSE; 481 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); 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) { 642 PetscCall(DMGetNamedGlobalVector(dm, "TSGLEE_Ydot", Ydot)); 643 } else *Ydot = glee->Ydot; 644 } 645 PetscFunctionReturn(PETSC_SUCCESS); 646 } 647 648 static PetscErrorCode TSGLEERestoreVecs(TS ts, DM dm, Vec *Ydot) 649 { 650 PetscFunctionBegin; 651 if (Ydot) { 652 if (dm && dm != ts->dm) PetscCall(DMRestoreNamedGlobalVector(dm, "TSGLEE_Ydot", Ydot)); 653 } 654 PetscFunctionReturn(PETSC_SUCCESS); 655 } 656 657 /* 658 This defines the nonlinear equation that is to be solved with SNES 659 */ 660 static PetscErrorCode SNESTSFormFunction_GLEE(SNES snes, Vec X, Vec F, TS ts) 661 { 662 TS_GLEE *glee = (TS_GLEE *)ts->data; 663 DM dm, dmsave; 664 Vec Ydot; 665 PetscReal shift = glee->scoeff / ts->time_step; 666 667 PetscFunctionBegin; 668 PetscCall(SNESGetDM(snes, &dm)); 669 PetscCall(TSGLEEGetVecs(ts, dm, &Ydot)); 670 /* Set Ydot = shift*X */ 671 PetscCall(VecCopy(X, Ydot)); 672 PetscCall(VecScale(Ydot, shift)); 673 dmsave = ts->dm; 674 ts->dm = dm; 675 676 PetscCall(TSComputeIFunction(ts, glee->stage_time, X, Ydot, F, PETSC_FALSE)); 677 678 ts->dm = dmsave; 679 PetscCall(TSGLEERestoreVecs(ts, dm, &Ydot)); 680 PetscFunctionReturn(PETSC_SUCCESS); 681 } 682 683 static PetscErrorCode SNESTSFormJacobian_GLEE(SNES snes, Vec X, Mat A, Mat B, TS ts) 684 { 685 TS_GLEE *glee = (TS_GLEE *)ts->data; 686 DM dm, dmsave; 687 Vec Ydot; 688 PetscReal shift = glee->scoeff / ts->time_step; 689 690 PetscFunctionBegin; 691 PetscCall(SNESGetDM(snes, &dm)); 692 PetscCall(TSGLEEGetVecs(ts, dm, &Ydot)); 693 /* glee->Ydot has already been computed in SNESTSFormFunction_GLEE (SNES guarantees this) */ 694 dmsave = ts->dm; 695 ts->dm = dm; 696 697 PetscCall(TSComputeIJacobian(ts, glee->stage_time, X, Ydot, shift, A, B, PETSC_FALSE)); 698 699 ts->dm = dmsave; 700 PetscCall(TSGLEERestoreVecs(ts, dm, &Ydot)); 701 PetscFunctionReturn(PETSC_SUCCESS); 702 } 703 704 static PetscErrorCode DMCoarsenHook_TSGLEE(DM fine, DM coarse, void *ctx) 705 { 706 PetscFunctionBegin; 707 PetscFunctionReturn(PETSC_SUCCESS); 708 } 709 710 static PetscErrorCode DMRestrictHook_TSGLEE(DM fine, Mat restrct, Vec rscale, Mat inject, DM coarse, void *ctx) 711 { 712 PetscFunctionBegin; 713 PetscFunctionReturn(PETSC_SUCCESS); 714 } 715 716 static PetscErrorCode DMSubDomainHook_TSGLEE(DM dm, DM subdm, void *ctx) 717 { 718 PetscFunctionBegin; 719 PetscFunctionReturn(PETSC_SUCCESS); 720 } 721 722 static PetscErrorCode DMSubDomainRestrictHook_TSGLEE(DM dm, VecScatter gscat, VecScatter lscat, DM subdm, void *ctx) 723 { 724 PetscFunctionBegin; 725 PetscFunctionReturn(PETSC_SUCCESS); 726 } 727 728 static PetscErrorCode TSSetUp_GLEE(TS ts) 729 { 730 TS_GLEE *glee = (TS_GLEE *)ts->data; 731 GLEETableau tab; 732 PetscInt s, r; 733 DM dm; 734 735 PetscFunctionBegin; 736 if (!glee->tableau) PetscCall(TSGLEESetType(ts, TSGLEEDefaultType)); 737 tab = glee->tableau; 738 s = tab->s; 739 r = tab->r; 740 PetscCall(VecDuplicateVecs(ts->vec_sol, r, &glee->Y)); 741 PetscCall(VecDuplicateVecs(ts->vec_sol, r, &glee->X)); 742 PetscCall(VecDuplicateVecs(ts->vec_sol, s, &glee->YStage)); 743 PetscCall(VecDuplicateVecs(ts->vec_sol, s, &glee->YdotStage)); 744 PetscCall(VecDuplicate(ts->vec_sol, &glee->Ydot)); 745 PetscCall(VecDuplicate(ts->vec_sol, &glee->yGErr)); 746 PetscCall(VecZeroEntries(glee->yGErr)); 747 PetscCall(VecDuplicate(ts->vec_sol, &glee->W)); 748 PetscCall(PetscMalloc2(s, &glee->swork, r, &glee->rwork)); 749 PetscCall(TSGetDM(ts, &dm)); 750 PetscCall(DMCoarsenHookAdd(dm, DMCoarsenHook_TSGLEE, DMRestrictHook_TSGLEE, ts)); 751 PetscCall(DMSubDomainHookAdd(dm, DMSubDomainHook_TSGLEE, DMSubDomainRestrictHook_TSGLEE, ts)); 752 PetscFunctionReturn(PETSC_SUCCESS); 753 } 754 755 static PetscErrorCode TSStartingMethod_GLEE(TS ts) 756 { 757 TS_GLEE *glee = (TS_GLEE *)ts->data; 758 GLEETableau tab = glee->tableau; 759 PetscInt r = tab->r, i; 760 PetscReal *S = tab->S; 761 762 PetscFunctionBegin; 763 for (i = 0; i < r; i++) { 764 PetscCall(VecZeroEntries(glee->Y[i])); 765 PetscCall(VecAXPY(glee->Y[i], S[i], ts->vec_sol)); 766 } 767 PetscFunctionReturn(PETSC_SUCCESS); 768 } 769 770 /*------------------------------------------------------------*/ 771 772 static PetscErrorCode TSSetFromOptions_GLEE(TS ts, PetscOptionItems *PetscOptionsObject) 773 { 774 char gleetype[256]; 775 776 PetscFunctionBegin; 777 PetscOptionsHeadBegin(PetscOptionsObject, "GLEE ODE solver options"); 778 { 779 GLEETableauLink link; 780 PetscInt count, choice; 781 PetscBool flg; 782 const char **namelist; 783 784 PetscCall(PetscStrncpy(gleetype, TSGLEEDefaultType, sizeof(gleetype))); 785 for (link = GLEETableauList, count = 0; link; link = link->next, count++); 786 PetscCall(PetscMalloc1(count, (char ***)&namelist)); 787 for (link = GLEETableauList, count = 0; link; link = link->next, count++) namelist[count] = link->tab.name; 788 PetscCall(PetscOptionsEList("-ts_glee_type", "Family of GLEE method", "TSGLEESetType", (const char *const *)namelist, count, gleetype, &choice, &flg)); 789 PetscCall(TSGLEESetType(ts, flg ? namelist[choice] : gleetype)); 790 PetscCall(PetscFree(namelist)); 791 } 792 PetscOptionsHeadEnd(); 793 PetscFunctionReturn(PETSC_SUCCESS); 794 } 795 796 static PetscErrorCode TSView_GLEE(TS ts, PetscViewer viewer) 797 { 798 TS_GLEE *glee = (TS_GLEE *)ts->data; 799 GLEETableau tab = glee->tableau; 800 PetscBool iascii; 801 802 PetscFunctionBegin; 803 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii)); 804 if (iascii) { 805 TSGLEEType gleetype; 806 char buf[512]; 807 PetscCall(TSGLEEGetType(ts, &gleetype)); 808 PetscCall(PetscViewerASCIIPrintf(viewer, " GLEE type %s\n", gleetype)); 809 PetscCall(PetscFormatRealArray(buf, sizeof(buf), "% 8.6f", tab->s, tab->c)); 810 PetscCall(PetscViewerASCIIPrintf(viewer, " Abscissa c = %s\n", buf)); 811 /* Note: print out r as well */ 812 } 813 PetscFunctionReturn(PETSC_SUCCESS); 814 } 815 816 static PetscErrorCode TSLoad_GLEE(TS ts, PetscViewer viewer) 817 { 818 SNES snes; 819 TSAdapt tsadapt; 820 821 PetscFunctionBegin; 822 PetscCall(TSGetAdapt(ts, &tsadapt)); 823 PetscCall(TSAdaptLoad(tsadapt, viewer)); 824 PetscCall(TSGetSNES(ts, &snes)); 825 PetscCall(SNESLoad(snes, viewer)); 826 /* function and Jacobian context for SNES when used with TS is always ts object */ 827 PetscCall(SNESSetFunction(snes, NULL, NULL, ts)); 828 PetscCall(SNESSetJacobian(snes, NULL, NULL, NULL, ts)); 829 PetscFunctionReturn(PETSC_SUCCESS); 830 } 831 832 /*@C 833 TSGLEESetType - Set the type of `TSGLEE` scheme 834 835 Logically Collective 836 837 Input Parameters: 838 + ts - timestepping context 839 - gleetype - type of `TSGLEE` scheme 840 841 Level: intermediate 842 843 .seealso: [](ch_ts), `TSGLEEGetType()`, `TSGLEE` 844 @*/ 845 PetscErrorCode TSGLEESetType(TS ts, TSGLEEType gleetype) 846 { 847 PetscFunctionBegin; 848 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 849 PetscAssertPointer(gleetype, 2); 850 PetscTryMethod(ts, "TSGLEESetType_C", (TS, TSGLEEType), (ts, gleetype)); 851 PetscFunctionReturn(PETSC_SUCCESS); 852 } 853 854 /*@C 855 TSGLEEGetType - Get the type of `TSGLEE` scheme 856 857 Logically Collective 858 859 Input Parameter: 860 . ts - timestepping context 861 862 Output Parameter: 863 . gleetype - type of `TSGLEE` scheme 864 865 Level: intermediate 866 867 .seealso: [](ch_ts), `TSGLEE`, `TSGLEESetType()` 868 @*/ 869 PetscErrorCode TSGLEEGetType(TS ts, TSGLEEType *gleetype) 870 { 871 PetscFunctionBegin; 872 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 873 PetscUseMethod(ts, "TSGLEEGetType_C", (TS, TSGLEEType *), (ts, gleetype)); 874 PetscFunctionReturn(PETSC_SUCCESS); 875 } 876 877 static PetscErrorCode TSGLEEGetType_GLEE(TS ts, TSGLEEType *gleetype) 878 { 879 TS_GLEE *glee = (TS_GLEE *)ts->data; 880 881 PetscFunctionBegin; 882 if (!glee->tableau) PetscCall(TSGLEESetType(ts, TSGLEEDefaultType)); 883 *gleetype = glee->tableau->name; 884 PetscFunctionReturn(PETSC_SUCCESS); 885 } 886 static PetscErrorCode TSGLEESetType_GLEE(TS ts, TSGLEEType gleetype) 887 { 888 TS_GLEE *glee = (TS_GLEE *)ts->data; 889 PetscBool match; 890 GLEETableauLink link; 891 892 PetscFunctionBegin; 893 if (glee->tableau) { 894 PetscCall(PetscStrcmp(glee->tableau->name, gleetype, &match)); 895 if (match) PetscFunctionReturn(PETSC_SUCCESS); 896 } 897 for (link = GLEETableauList; link; link = link->next) { 898 PetscCall(PetscStrcmp(link->tab.name, gleetype, &match)); 899 if (match) { 900 PetscCall(TSReset_GLEE(ts)); 901 glee->tableau = &link->tab; 902 PetscFunctionReturn(PETSC_SUCCESS); 903 } 904 } 905 SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_UNKNOWN_TYPE, "Could not find '%s'", gleetype); 906 } 907 908 static PetscErrorCode TSGetStages_GLEE(TS ts, PetscInt *ns, Vec **Y) 909 { 910 TS_GLEE *glee = (TS_GLEE *)ts->data; 911 912 PetscFunctionBegin; 913 if (ns) *ns = glee->tableau->s; 914 if (Y) *Y = glee->YStage; 915 PetscFunctionReturn(PETSC_SUCCESS); 916 } 917 918 static PetscErrorCode TSGetSolutionComponents_GLEE(TS ts, PetscInt *n, Vec *Y) 919 { 920 TS_GLEE *glee = (TS_GLEE *)ts->data; 921 GLEETableau tab = glee->tableau; 922 923 PetscFunctionBegin; 924 if (!Y) *n = tab->r; 925 else { 926 if ((*n >= 0) && (*n < tab->r)) { 927 PetscCall(VecCopy(glee->Y[*n], *Y)); 928 } else SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_OUTOFRANGE, "Second argument (%" PetscInt_FMT ") out of range[0,%" PetscInt_FMT "].", *n, tab->r - 1); 929 } 930 PetscFunctionReturn(PETSC_SUCCESS); 931 } 932 933 static PetscErrorCode TSGetAuxSolution_GLEE(TS ts, Vec *X) 934 { 935 TS_GLEE *glee = (TS_GLEE *)ts->data; 936 GLEETableau tab = glee->tableau; 937 PetscReal *F = tab->Fembed; 938 PetscInt r = tab->r; 939 Vec *Y = glee->Y; 940 PetscScalar *wr = glee->rwork; 941 PetscInt i; 942 943 PetscFunctionBegin; 944 PetscCall(VecZeroEntries(*X)); 945 for (i = 0; i < r; i++) wr[i] = F[i]; 946 PetscCall(VecMAXPY(*X, r, wr, Y)); 947 PetscFunctionReturn(PETSC_SUCCESS); 948 } 949 950 static PetscErrorCode TSGetTimeError_GLEE(TS ts, PetscInt n, Vec *X) 951 { 952 TS_GLEE *glee = (TS_GLEE *)ts->data; 953 GLEETableau tab = glee->tableau; 954 PetscReal *F = tab->Ferror; 955 PetscInt r = tab->r; 956 Vec *Y = glee->Y; 957 PetscScalar *wr = glee->rwork; 958 PetscInt i; 959 960 PetscFunctionBegin; 961 PetscCall(VecZeroEntries(*X)); 962 if (n == 0) { 963 for (i = 0; i < r; i++) wr[i] = F[i]; 964 PetscCall(VecMAXPY(*X, r, wr, Y)); 965 } else if (n == -1) { 966 *X = glee->yGErr; 967 } 968 PetscFunctionReturn(PETSC_SUCCESS); 969 } 970 971 static PetscErrorCode TSSetTimeError_GLEE(TS ts, Vec X) 972 { 973 TS_GLEE *glee = (TS_GLEE *)ts->data; 974 GLEETableau tab = glee->tableau; 975 PetscReal *S = tab->Serror; 976 PetscInt r = tab->r, i; 977 Vec *Y = glee->Y; 978 979 PetscFunctionBegin; 980 PetscCheck(r == 2, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "TSSetTimeError_GLEE not supported for '%s' with r=%" PetscInt_FMT ".", tab->name, tab->r); 981 for (i = 1; i < r; i++) { 982 PetscCall(VecCopy(ts->vec_sol, Y[i])); 983 PetscCall(VecAXPBY(Y[i], S[0], S[1], X)); 984 PetscCall(VecCopy(X, glee->yGErr)); 985 } 986 PetscFunctionReturn(PETSC_SUCCESS); 987 } 988 989 static PetscErrorCode TSDestroy_GLEE(TS ts) 990 { 991 PetscFunctionBegin; 992 PetscCall(TSReset_GLEE(ts)); 993 if (ts->dm) { 994 PetscCall(DMCoarsenHookRemove(ts->dm, DMCoarsenHook_TSGLEE, DMRestrictHook_TSGLEE, ts)); 995 PetscCall(DMSubDomainHookRemove(ts->dm, DMSubDomainHook_TSGLEE, DMSubDomainRestrictHook_TSGLEE, ts)); 996 } 997 PetscCall(PetscFree(ts->data)); 998 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSGLEEGetType_C", NULL)); 999 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSGLEESetType_C", NULL)); 1000 PetscFunctionReturn(PETSC_SUCCESS); 1001 } 1002 1003 /* ------------------------------------------------------------ */ 1004 /*MC 1005 TSGLEE - ODE and DAE solver using General Linear with Error Estimation schemes 1006 1007 The user should provide the right-hand side of the equation using `TSSetRHSFunction()`. 1008 1009 Level: beginner 1010 1011 Note: 1012 The default is `TSGLEE35`, it can be changed with `TSGLEESetType()` or -ts_glee_type 1013 1014 .seealso: [](ch_ts), `TSCreate()`, `TS`, `TSSetType()`, `TSGLEESetType()`, `TSGLEEGetType()`, 1015 `TSGLEE23`, `TTSGLEE24`, `TSGLEE35`, `TSGLEE25I`, `TSGLEEEXRK2A`, 1016 `TSGLEERK32G1`, `TSGLEERK285EX`, `TSGLEEType`, `TSGLEERegister()`, `TSType` 1017 M*/ 1018 PETSC_EXTERN PetscErrorCode TSCreate_GLEE(TS ts) 1019 { 1020 TS_GLEE *th; 1021 1022 PetscFunctionBegin; 1023 PetscCall(TSGLEEInitializePackage()); 1024 1025 ts->ops->reset = TSReset_GLEE; 1026 ts->ops->destroy = TSDestroy_GLEE; 1027 ts->ops->view = TSView_GLEE; 1028 ts->ops->load = TSLoad_GLEE; 1029 ts->ops->setup = TSSetUp_GLEE; 1030 ts->ops->step = TSStep_GLEE; 1031 ts->ops->interpolate = TSInterpolate_GLEE; 1032 ts->ops->evaluatestep = TSEvaluateStep_GLEE; 1033 ts->ops->setfromoptions = TSSetFromOptions_GLEE; 1034 ts->ops->getstages = TSGetStages_GLEE; 1035 ts->ops->snesfunction = SNESTSFormFunction_GLEE; 1036 ts->ops->snesjacobian = SNESTSFormJacobian_GLEE; 1037 ts->ops->getsolutioncomponents = TSGetSolutionComponents_GLEE; 1038 ts->ops->getauxsolution = TSGetAuxSolution_GLEE; 1039 ts->ops->gettimeerror = TSGetTimeError_GLEE; 1040 ts->ops->settimeerror = TSSetTimeError_GLEE; 1041 ts->ops->startingmethod = TSStartingMethod_GLEE; 1042 ts->default_adapt_type = TSADAPTGLEE; 1043 1044 ts->usessnes = PETSC_TRUE; 1045 1046 PetscCall(PetscNew(&th)); 1047 ts->data = (void *)th; 1048 1049 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSGLEEGetType_C", TSGLEEGetType_GLEE)); 1050 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSGLEESetType_C", TSGLEESetType_GLEE)); 1051 PetscFunctionReturn(PETSC_SUCCESS); 1052 } 1053