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