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 ; 787 PetscCall(PetscMalloc1(count, (char ***)&namelist)); 788 for (link = GLEETableauList, count = 0; link; link = link->next, count++) namelist[count] = link->tab.name; 789 PetscCall(PetscOptionsEList("-ts_glee_type", "Family of GLEE method", "TSGLEESetType", (const char *const *)namelist, count, gleetype, &choice, &flg)); 790 PetscCall(TSGLEESetType(ts, flg ? namelist[choice] : gleetype)); 791 PetscCall(PetscFree(namelist)); 792 } 793 PetscOptionsHeadEnd(); 794 PetscFunctionReturn(PETSC_SUCCESS); 795 } 796 797 static PetscErrorCode TSView_GLEE(TS ts, PetscViewer viewer) 798 { 799 TS_GLEE *glee = (TS_GLEE *)ts->data; 800 GLEETableau tab = glee->tableau; 801 PetscBool iascii; 802 803 PetscFunctionBegin; 804 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii)); 805 if (iascii) { 806 TSGLEEType gleetype; 807 char buf[512]; 808 PetscCall(TSGLEEGetType(ts, &gleetype)); 809 PetscCall(PetscViewerASCIIPrintf(viewer, " GLEE type %s\n", gleetype)); 810 PetscCall(PetscFormatRealArray(buf, sizeof(buf), "% 8.6f", tab->s, tab->c)); 811 PetscCall(PetscViewerASCIIPrintf(viewer, " Abscissa c = %s\n", buf)); 812 /* Note: print out r as well */ 813 } 814 PetscFunctionReturn(PETSC_SUCCESS); 815 } 816 817 static PetscErrorCode TSLoad_GLEE(TS ts, PetscViewer viewer) 818 { 819 SNES snes; 820 TSAdapt tsadapt; 821 822 PetscFunctionBegin; 823 PetscCall(TSGetAdapt(ts, &tsadapt)); 824 PetscCall(TSAdaptLoad(tsadapt, viewer)); 825 PetscCall(TSGetSNES(ts, &snes)); 826 PetscCall(SNESLoad(snes, viewer)); 827 /* function and Jacobian context for SNES when used with TS is always ts object */ 828 PetscCall(SNESSetFunction(snes, NULL, NULL, ts)); 829 PetscCall(SNESSetJacobian(snes, NULL, NULL, NULL, ts)); 830 PetscFunctionReturn(PETSC_SUCCESS); 831 } 832 833 /*@C 834 TSGLEESetType - Set the type of `TSGLEE` scheme 835 836 Logically Collective 837 838 Input Parameters: 839 + ts - timestepping context 840 - gleetype - type of `TSGLEE` scheme 841 842 Level: intermediate 843 844 .seealso: [](ch_ts), `TSGLEEGetType()`, `TSGLEE` 845 @*/ 846 PetscErrorCode TSGLEESetType(TS ts, TSGLEEType gleetype) 847 { 848 PetscFunctionBegin; 849 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 850 PetscAssertPointer(gleetype, 2); 851 PetscTryMethod(ts, "TSGLEESetType_C", (TS, TSGLEEType), (ts, gleetype)); 852 PetscFunctionReturn(PETSC_SUCCESS); 853 } 854 855 /*@C 856 TSGLEEGetType - Get the type of `TSGLEE` scheme 857 858 Logically Collective 859 860 Input Parameter: 861 . ts - timestepping context 862 863 Output Parameter: 864 . gleetype - type of `TSGLEE` scheme 865 866 Level: intermediate 867 868 .seealso: [](ch_ts), `TSGLEE`, `TSGLEESetType()` 869 @*/ 870 PetscErrorCode TSGLEEGetType(TS ts, TSGLEEType *gleetype) 871 { 872 PetscFunctionBegin; 873 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 874 PetscUseMethod(ts, "TSGLEEGetType_C", (TS, TSGLEEType *), (ts, gleetype)); 875 PetscFunctionReturn(PETSC_SUCCESS); 876 } 877 878 static PetscErrorCode TSGLEEGetType_GLEE(TS ts, TSGLEEType *gleetype) 879 { 880 TS_GLEE *glee = (TS_GLEE *)ts->data; 881 882 PetscFunctionBegin; 883 if (!glee->tableau) PetscCall(TSGLEESetType(ts, TSGLEEDefaultType)); 884 *gleetype = glee->tableau->name; 885 PetscFunctionReturn(PETSC_SUCCESS); 886 } 887 static PetscErrorCode TSGLEESetType_GLEE(TS ts, TSGLEEType gleetype) 888 { 889 TS_GLEE *glee = (TS_GLEE *)ts->data; 890 PetscBool match; 891 GLEETableauLink link; 892 893 PetscFunctionBegin; 894 if (glee->tableau) { 895 PetscCall(PetscStrcmp(glee->tableau->name, gleetype, &match)); 896 if (match) PetscFunctionReturn(PETSC_SUCCESS); 897 } 898 for (link = GLEETableauList; link; link = link->next) { 899 PetscCall(PetscStrcmp(link->tab.name, gleetype, &match)); 900 if (match) { 901 PetscCall(TSReset_GLEE(ts)); 902 glee->tableau = &link->tab; 903 PetscFunctionReturn(PETSC_SUCCESS); 904 } 905 } 906 SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_UNKNOWN_TYPE, "Could not find '%s'", gleetype); 907 } 908 909 static PetscErrorCode TSGetStages_GLEE(TS ts, PetscInt *ns, Vec **Y) 910 { 911 TS_GLEE *glee = (TS_GLEE *)ts->data; 912 913 PetscFunctionBegin; 914 if (ns) *ns = glee->tableau->s; 915 if (Y) *Y = glee->YStage; 916 PetscFunctionReturn(PETSC_SUCCESS); 917 } 918 919 static PetscErrorCode TSGetSolutionComponents_GLEE(TS ts, PetscInt *n, Vec *Y) 920 { 921 TS_GLEE *glee = (TS_GLEE *)ts->data; 922 GLEETableau tab = glee->tableau; 923 924 PetscFunctionBegin; 925 if (!Y) *n = tab->r; 926 else { 927 if ((*n >= 0) && (*n < tab->r)) { 928 PetscCall(VecCopy(glee->Y[*n], *Y)); 929 } else SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_OUTOFRANGE, "Second argument (%" PetscInt_FMT ") out of range[0,%" PetscInt_FMT "].", *n, tab->r - 1); 930 } 931 PetscFunctionReturn(PETSC_SUCCESS); 932 } 933 934 static PetscErrorCode TSGetAuxSolution_GLEE(TS ts, Vec *X) 935 { 936 TS_GLEE *glee = (TS_GLEE *)ts->data; 937 GLEETableau tab = glee->tableau; 938 PetscReal *F = tab->Fembed; 939 PetscInt r = tab->r; 940 Vec *Y = glee->Y; 941 PetscScalar *wr = glee->rwork; 942 PetscInt i; 943 944 PetscFunctionBegin; 945 PetscCall(VecZeroEntries(*X)); 946 for (i = 0; i < r; i++) wr[i] = F[i]; 947 PetscCall(VecMAXPY(*X, r, wr, Y)); 948 PetscFunctionReturn(PETSC_SUCCESS); 949 } 950 951 static PetscErrorCode TSGetTimeError_GLEE(TS ts, PetscInt n, Vec *X) 952 { 953 TS_GLEE *glee = (TS_GLEE *)ts->data; 954 GLEETableau tab = glee->tableau; 955 PetscReal *F = tab->Ferror; 956 PetscInt r = tab->r; 957 Vec *Y = glee->Y; 958 PetscScalar *wr = glee->rwork; 959 PetscInt i; 960 961 PetscFunctionBegin; 962 PetscCall(VecZeroEntries(*X)); 963 if (n == 0) { 964 for (i = 0; i < r; i++) wr[i] = F[i]; 965 PetscCall(VecMAXPY(*X, r, wr, Y)); 966 } else if (n == -1) { 967 *X = glee->yGErr; 968 } 969 PetscFunctionReturn(PETSC_SUCCESS); 970 } 971 972 static PetscErrorCode TSSetTimeError_GLEE(TS ts, Vec X) 973 { 974 TS_GLEE *glee = (TS_GLEE *)ts->data; 975 GLEETableau tab = glee->tableau; 976 PetscReal *S = tab->Serror; 977 PetscInt r = tab->r, i; 978 Vec *Y = glee->Y; 979 980 PetscFunctionBegin; 981 PetscCheck(r == 2, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "TSSetTimeError_GLEE not supported for '%s' with r=%" PetscInt_FMT ".", tab->name, tab->r); 982 for (i = 1; i < r; i++) { 983 PetscCall(VecCopy(ts->vec_sol, Y[i])); 984 PetscCall(VecAXPBY(Y[i], S[0], S[1], X)); 985 PetscCall(VecCopy(X, glee->yGErr)); 986 } 987 PetscFunctionReturn(PETSC_SUCCESS); 988 } 989 990 static PetscErrorCode TSDestroy_GLEE(TS ts) 991 { 992 PetscFunctionBegin; 993 PetscCall(TSReset_GLEE(ts)); 994 if (ts->dm) { 995 PetscCall(DMCoarsenHookRemove(ts->dm, DMCoarsenHook_TSGLEE, DMRestrictHook_TSGLEE, ts)); 996 PetscCall(DMSubDomainHookRemove(ts->dm, DMSubDomainHook_TSGLEE, DMSubDomainRestrictHook_TSGLEE, ts)); 997 } 998 PetscCall(PetscFree(ts->data)); 999 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSGLEEGetType_C", NULL)); 1000 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSGLEESetType_C", NULL)); 1001 PetscFunctionReturn(PETSC_SUCCESS); 1002 } 1003 1004 /* ------------------------------------------------------------ */ 1005 /*MC 1006 TSGLEE - ODE and DAE solver using General Linear with Error Estimation schemes 1007 1008 The user should provide the right hand side of the equation using `TSSetRHSFunction()`. 1009 1010 Level: beginner 1011 1012 Note: 1013 The default is `TSGLEE35`, it can be changed with `TSGLEESetType()` or -ts_glee_type 1014 1015 .seealso: [](ch_ts), `TSCreate()`, `TS`, `TSSetType()`, `TSGLEESetType()`, `TSGLEEGetType()`, 1016 `TSGLEE23`, `TTSGLEE24`, `TSGLEE35`, `TSGLEE25I`, `TSGLEEEXRK2A`, 1017 `TSGLEERK32G1`, `TSGLEERK285EX`, `TSGLEEType`, `TSGLEERegister()`, `TSType` 1018 M*/ 1019 PETSC_EXTERN PetscErrorCode TSCreate_GLEE(TS ts) 1020 { 1021 TS_GLEE *th; 1022 1023 PetscFunctionBegin; 1024 PetscCall(TSGLEEInitializePackage()); 1025 1026 ts->ops->reset = TSReset_GLEE; 1027 ts->ops->destroy = TSDestroy_GLEE; 1028 ts->ops->view = TSView_GLEE; 1029 ts->ops->load = TSLoad_GLEE; 1030 ts->ops->setup = TSSetUp_GLEE; 1031 ts->ops->step = TSStep_GLEE; 1032 ts->ops->interpolate = TSInterpolate_GLEE; 1033 ts->ops->evaluatestep = TSEvaluateStep_GLEE; 1034 ts->ops->setfromoptions = TSSetFromOptions_GLEE; 1035 ts->ops->getstages = TSGetStages_GLEE; 1036 ts->ops->snesfunction = SNESTSFormFunction_GLEE; 1037 ts->ops->snesjacobian = SNESTSFormJacobian_GLEE; 1038 ts->ops->getsolutioncomponents = TSGetSolutionComponents_GLEE; 1039 ts->ops->getauxsolution = TSGetAuxSolution_GLEE; 1040 ts->ops->gettimeerror = TSGetTimeError_GLEE; 1041 ts->ops->settimeerror = TSSetTimeError_GLEE; 1042 ts->ops->startingmethod = TSStartingMethod_GLEE; 1043 ts->default_adapt_type = TSADAPTGLEE; 1044 1045 ts->usessnes = PETSC_TRUE; 1046 1047 PetscCall(PetscNew(&th)); 1048 ts->data = (void *)th; 1049 1050 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSGLEEGetType_C", TSGLEEGetType_GLEE)); 1051 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSGLEESetType_C", TSGLEESetType_GLEE)); 1052 PetscFunctionReturn(PETSC_SUCCESS); 1053 } 1054