1 #include <../src/ts/impls/implicit/glle/glle.h> /*I "petscts.h" I*/ 2 #include <petscdm.h> 3 #include <petscblaslapack.h> 4 5 static const char *TSGLLEErrorDirections[] = {"FORWARD", "BACKWARD", "TSGLLEErrorDirection", "TSGLLEERROR_", NULL}; 6 static PetscFunctionList TSGLLEList; 7 static PetscFunctionList TSGLLEAcceptList; 8 static PetscBool TSGLLEPackageInitialized; 9 static PetscBool TSGLLERegisterAllCalled; 10 11 /* This function is pure */ 12 static PetscScalar Factorial(PetscInt n) 13 { 14 PetscInt i; 15 if (n < 12) { /* Can compute with 32-bit integers */ 16 PetscInt f = 1; 17 for (i = 2; i <= n; i++) f *= i; 18 return (PetscScalar)f; 19 } else { 20 PetscScalar f = 1.; 21 for (i = 2; i <= n; i++) f *= (PetscScalar)i; 22 return f; 23 } 24 } 25 26 /* This function is pure */ 27 static PetscScalar CPowF(PetscScalar c, PetscInt p) 28 { 29 return PetscPowRealInt(PetscRealPart(c), p) / Factorial(p); 30 } 31 32 static PetscErrorCode TSGLLEGetVecs(TS ts, DM dm, Vec *Z, Vec *Ydotstage) 33 { 34 TS_GLLE *gl = (TS_GLLE *)ts->data; 35 36 PetscFunctionBegin; 37 if (Z) { 38 if (dm && dm != ts->dm) { 39 PetscCall(DMGetNamedGlobalVector(dm, "TSGLLE_Z", Z)); 40 } else *Z = gl->Z; 41 } 42 if (Ydotstage) { 43 if (dm && dm != ts->dm) { 44 PetscCall(DMGetNamedGlobalVector(dm, "TSGLLE_Ydot", Ydotstage)); 45 } else *Ydotstage = gl->Ydot[gl->stage]; 46 } 47 PetscFunctionReturn(PETSC_SUCCESS); 48 } 49 50 static PetscErrorCode TSGLLERestoreVecs(TS ts, DM dm, Vec *Z, Vec *Ydotstage) 51 { 52 PetscFunctionBegin; 53 if (Z) { 54 if (dm && dm != ts->dm) PetscCall(DMRestoreNamedGlobalVector(dm, "TSGLLE_Z", Z)); 55 } 56 if (Ydotstage) { 57 if (dm && dm != ts->dm) PetscCall(DMRestoreNamedGlobalVector(dm, "TSGLLE_Ydot", Ydotstage)); 58 } 59 PetscFunctionReturn(PETSC_SUCCESS); 60 } 61 62 static PetscErrorCode DMCoarsenHook_TSGLLE(DM fine, DM coarse, void *ctx) 63 { 64 PetscFunctionBegin; 65 PetscFunctionReturn(PETSC_SUCCESS); 66 } 67 68 static PetscErrorCode DMRestrictHook_TSGLLE(DM fine, Mat restrct, Vec rscale, Mat inject, DM coarse, void *ctx) 69 { 70 TS ts = (TS)ctx; 71 Vec Ydot, Ydot_c; 72 73 PetscFunctionBegin; 74 PetscCall(TSGLLEGetVecs(ts, fine, NULL, &Ydot)); 75 PetscCall(TSGLLEGetVecs(ts, coarse, NULL, &Ydot_c)); 76 PetscCall(MatRestrict(restrct, Ydot, Ydot_c)); 77 PetscCall(VecPointwiseMult(Ydot_c, rscale, Ydot_c)); 78 PetscCall(TSGLLERestoreVecs(ts, fine, NULL, &Ydot)); 79 PetscCall(TSGLLERestoreVecs(ts, coarse, NULL, &Ydot_c)); 80 PetscFunctionReturn(PETSC_SUCCESS); 81 } 82 83 static PetscErrorCode DMSubDomainHook_TSGLLE(DM dm, DM subdm, void *ctx) 84 { 85 PetscFunctionBegin; 86 PetscFunctionReturn(PETSC_SUCCESS); 87 } 88 89 static PetscErrorCode DMSubDomainRestrictHook_TSGLLE(DM dm, VecScatter gscat, VecScatter lscat, DM subdm, void *ctx) 90 { 91 TS ts = (TS)ctx; 92 Vec Ydot, Ydot_s; 93 94 PetscFunctionBegin; 95 PetscCall(TSGLLEGetVecs(ts, dm, NULL, &Ydot)); 96 PetscCall(TSGLLEGetVecs(ts, subdm, NULL, &Ydot_s)); 97 98 PetscCall(VecScatterBegin(gscat, Ydot, Ydot_s, INSERT_VALUES, SCATTER_FORWARD)); 99 PetscCall(VecScatterEnd(gscat, Ydot, Ydot_s, INSERT_VALUES, SCATTER_FORWARD)); 100 101 PetscCall(TSGLLERestoreVecs(ts, dm, NULL, &Ydot)); 102 PetscCall(TSGLLERestoreVecs(ts, subdm, NULL, &Ydot_s)); 103 PetscFunctionReturn(PETSC_SUCCESS); 104 } 105 106 static PetscErrorCode TSGLLESchemeCreate(PetscInt p, PetscInt q, PetscInt r, PetscInt s, const PetscScalar *c, const PetscScalar *a, const PetscScalar *b, const PetscScalar *u, const PetscScalar *v, TSGLLEScheme *inscheme) 107 { 108 TSGLLEScheme scheme; 109 PetscInt j; 110 111 PetscFunctionBegin; 112 PetscCheck(p >= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Scheme order must be positive"); 113 PetscCheck(r >= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "At least one item must be carried between steps"); 114 PetscCheck(s >= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "At least one stage is required"); 115 PetscAssertPointer(inscheme, 10); 116 *inscheme = NULL; 117 PetscCall(PetscNew(&scheme)); 118 scheme->p = p; 119 scheme->q = q; 120 scheme->r = r; 121 scheme->s = s; 122 123 PetscCall(PetscMalloc5(s, &scheme->c, s * s, &scheme->a, r * s, &scheme->b, r * s, &scheme->u, r * r, &scheme->v)); 124 PetscCall(PetscArraycpy(scheme->c, c, s)); 125 for (j = 0; j < s * s; j++) scheme->a[j] = (PetscAbsScalar(a[j]) < 1e-12) ? 0 : a[j]; 126 for (j = 0; j < r * s; j++) scheme->b[j] = (PetscAbsScalar(b[j]) < 1e-12) ? 0 : b[j]; 127 for (j = 0; j < s * r; j++) scheme->u[j] = (PetscAbsScalar(u[j]) < 1e-12) ? 0 : u[j]; 128 for (j = 0; j < r * r; j++) scheme->v[j] = (PetscAbsScalar(v[j]) < 1e-12) ? 0 : v[j]; 129 130 PetscCall(PetscMalloc6(r, &scheme->alpha, r, &scheme->beta, r, &scheme->gamma, 3 * s, &scheme->phi, 3 * r, &scheme->psi, r, &scheme->stage_error)); 131 { 132 PetscInt i, j, k, ss = s + 2; 133 PetscBLASInt m, n, one = 1, *ipiv, lwork, info, ldb; 134 PetscReal rcond, *sing, *workreal; 135 PetscScalar *ImV, *H, *bmat, *workscalar, *c = scheme->c, *a = scheme->a, *b = scheme->b, *u = scheme->u, *v = scheme->v; 136 PetscBLASInt rank; 137 138 PetscCall(PetscBLASIntCast(4 * ((s + 3) * 3 + 3), &lwork)); 139 PetscCall(PetscMalloc7(PetscSqr(r), &ImV, 3 * s, &H, 3 * ss, &bmat, lwork, &workscalar, 5 * (3 + r), &workreal, r + s, &sing, r + s, &ipiv)); 140 141 /* column-major input */ 142 for (i = 0; i < r - 1; i++) { 143 for (j = 0; j < r - 1; j++) ImV[i + j * r] = 1.0 * (i == j) - v[(i + 1) * r + j + 1]; 144 } 145 /* Build right-hand side for alpha (tp - glm.B(2:end,:)*(glm.c.^(p)./factorial(p))) */ 146 for (i = 1; i < r; i++) { 147 scheme->alpha[i] = 1. / Factorial(p + 1 - i); 148 for (j = 0; j < s; j++) scheme->alpha[i] -= b[i * s + j] * CPowF(c[j], p); 149 } 150 PetscCall(PetscBLASIntCast(r - 1, &m)); 151 PetscCall(PetscBLASIntCast(r, &n)); 152 PetscCallBLAS("LAPACKgesv", LAPACKgesv_(&m, &one, ImV, &n, ipiv, scheme->alpha + 1, &n, &info)); 153 PetscCheck(info >= 0, PETSC_COMM_SELF, PETSC_ERR_LIB, "Bad argument to GESV"); 154 PetscCheck(info <= 0, PETSC_COMM_SELF, PETSC_ERR_MAT_LU_ZRPVT, "Bad LU factorization"); 155 156 /* Build right-hand side for beta (tp1 - glm.B(2:end,:)*(glm.c.^(p+1)./factorial(p+1)) - e.alpha) */ 157 for (i = 1; i < r; i++) { 158 scheme->beta[i] = 1. / Factorial(p + 2 - i) - scheme->alpha[i]; 159 for (j = 0; j < s; j++) scheme->beta[i] -= b[i * s + j] * CPowF(c[j], p + 1); 160 } 161 PetscCallBLAS("LAPACKgetrs", LAPACKgetrs_("No transpose", &m, &one, ImV, &n, ipiv, scheme->beta + 1, &n, &info)); 162 PetscCheck(info >= 0, PETSC_COMM_SELF, PETSC_ERR_LIB, "Bad argument to GETRS"); 163 PetscCheck(info <= 0, PETSC_COMM_SELF, PETSC_ERR_LIB, "Should not happen"); 164 165 /* Build stage_error vector 166 xi = glm.c.^(p+1)/factorial(p+1) - glm.A*glm.c.^p/factorial(p) + glm.U(:,2:end)*e.alpha; 167 */ 168 for (i = 0; i < s; i++) { 169 scheme->stage_error[i] = CPowF(c[i], p + 1); 170 for (j = 0; j < s; j++) scheme->stage_error[i] -= a[i * s + j] * CPowF(c[j], p); 171 for (j = 1; j < r; j++) scheme->stage_error[i] += u[i * r + j] * scheme->alpha[j]; 172 } 173 174 /* alpha[0] (epsilon in B,J,W 2007) 175 epsilon = 1/factorial(p+1) - B(1,:)*c.^p/factorial(p) + V(1,2:end)*e.alpha; 176 */ 177 scheme->alpha[0] = 1. / Factorial(p + 1); 178 for (j = 0; j < s; j++) scheme->alpha[0] -= b[0 * s + j] * CPowF(c[j], p); 179 for (j = 1; j < r; j++) scheme->alpha[0] += v[0 * r + j] * scheme->alpha[j]; 180 181 /* right-hand side for gamma (glm.B(2:end,:)*e.xi - e.epsilon*eye(s-1,1)) */ 182 for (i = 1; i < r; i++) { 183 scheme->gamma[i] = (i == 1 ? -1. : 0) * scheme->alpha[0]; 184 for (j = 0; j < s; j++) scheme->gamma[i] += b[i * s + j] * scheme->stage_error[j]; 185 } 186 PetscCallBLAS("LAPACKgetrs", LAPACKgetrs_("No transpose", &m, &one, ImV, &n, ipiv, scheme->gamma + 1, &n, &info)); 187 PetscCheck(info >= 0, PETSC_COMM_SELF, PETSC_ERR_LIB, "Bad argument to GETRS"); 188 PetscCheck(info <= 0, PETSC_COMM_SELF, PETSC_ERR_LIB, "Should not happen"); 189 190 /* beta[0] (rho in B,J,W 2007) 191 e.rho = 1/factorial(p+2) - glm.B(1,:)*glm.c.^(p+1)/factorial(p+1) ... 192 + glm.V(1,2:end)*e.beta;% - e.epsilon; 193 % Note: The paper (B,J,W 2007) includes the last term in their definition 194 * */ 195 scheme->beta[0] = 1. / Factorial(p + 2); 196 for (j = 0; j < s; j++) scheme->beta[0] -= b[0 * s + j] * CPowF(c[j], p + 1); 197 for (j = 1; j < r; j++) scheme->beta[0] += v[0 * r + j] * scheme->beta[j]; 198 199 /* gamma[0] (sigma in B,J,W 2007) 200 * e.sigma = glm.B(1,:)*e.xi + glm.V(1,2:end)*e.gamma; 201 * */ 202 scheme->gamma[0] = 0.0; 203 for (j = 0; j < s; j++) scheme->gamma[0] += b[0 * s + j] * scheme->stage_error[j]; 204 for (j = 1; j < r; j++) scheme->gamma[0] += v[0 * s + j] * scheme->gamma[j]; 205 206 /* Assemble H 207 * % " PetscInt_FMT "etermine the error estimators phi 208 H = [[cpow(glm.c,p) + C*e.alpha] [cpow(glm.c,p+1) + C*e.beta] ... 209 [e.xi - C*(e.gamma + 0*e.epsilon*eye(s-1,1))]]'; 210 % Paper has formula above without the 0, but that term must be left 211 % out to satisfy the conditions they propose and to make the 212 % example schemes work 213 e.H = H; 214 e.phi = (H \ [1 0 0;1 1 0;0 0 -1])'; 215 e.psi = -e.phi*C; 216 * */ 217 for (j = 0; j < s; j++) { 218 H[0 + j * 3] = CPowF(c[j], p); 219 H[1 + j * 3] = CPowF(c[j], p + 1); 220 H[2 + j * 3] = scheme->stage_error[j]; 221 for (k = 1; k < r; k++) { 222 H[0 + j * 3] += CPowF(c[j], k - 1) * scheme->alpha[k]; 223 H[1 + j * 3] += CPowF(c[j], k - 1) * scheme->beta[k]; 224 H[2 + j * 3] -= CPowF(c[j], k - 1) * scheme->gamma[k]; 225 } 226 } 227 bmat[0 + 0 * ss] = 1.; 228 bmat[0 + 1 * ss] = 0.; 229 bmat[0 + 2 * ss] = 0.; 230 bmat[1 + 0 * ss] = 1.; 231 bmat[1 + 1 * ss] = 1.; 232 bmat[1 + 2 * ss] = 0.; 233 bmat[2 + 0 * ss] = 0.; 234 bmat[2 + 1 * ss] = 0.; 235 bmat[2 + 2 * ss] = -1.; 236 m = 3; 237 PetscCall(PetscBLASIntCast(s, &n)); 238 PetscCall(PetscBLASIntCast(ss, &ldb)); 239 rcond = 1e-12; 240 #if defined(PETSC_USE_COMPLEX) 241 /* ZGELSS( M, N, NRHS, A, LDA, B, LDB, S, RCOND, RANK, WORK, LWORK, RWORK, INFO) */ 242 PetscCallBLAS("LAPACKgelss", LAPACKgelss_(&m, &n, &m, H, &m, bmat, &ldb, sing, &rcond, &rank, workscalar, &lwork, workreal, &info)); 243 #else 244 /* DGELSS( M, N, NRHS, A, LDA, B, LDB, S, RCOND, RANK, WORK, LWORK, INFO) */ 245 PetscCallBLAS("LAPACKgelss", LAPACKgelss_(&m, &n, &m, H, &m, bmat, &ldb, sing, &rcond, &rank, workscalar, &lwork, &info)); 246 #endif 247 PetscCheck(info >= 0, PETSC_COMM_SELF, PETSC_ERR_LIB, "Bad argument to GELSS"); 248 PetscCheck(info <= 0, PETSC_COMM_SELF, PETSC_ERR_LIB, "SVD failed to converge"); 249 250 for (j = 0; j < 3; j++) { 251 for (k = 0; k < s; k++) scheme->phi[k + j * s] = bmat[k + j * ss]; 252 } 253 254 /* the other part of the error estimator, psi in B,J,W 2007 */ 255 scheme->psi[0 * r + 0] = 0.; 256 scheme->psi[1 * r + 0] = 0.; 257 scheme->psi[2 * r + 0] = 0.; 258 for (j = 1; j < r; j++) { 259 scheme->psi[0 * r + j] = 0.; 260 scheme->psi[1 * r + j] = 0.; 261 scheme->psi[2 * r + j] = 0.; 262 for (k = 0; k < s; k++) { 263 scheme->psi[0 * r + j] -= CPowF(c[k], j - 1) * scheme->phi[0 * s + k]; 264 scheme->psi[1 * r + j] -= CPowF(c[k], j - 1) * scheme->phi[1 * s + k]; 265 scheme->psi[2 * r + j] -= CPowF(c[k], j - 1) * scheme->phi[2 * s + k]; 266 } 267 } 268 PetscCall(PetscFree7(ImV, H, bmat, workscalar, workreal, sing, ipiv)); 269 } 270 /* Check which properties are satisfied */ 271 scheme->stiffly_accurate = PETSC_TRUE; 272 if (scheme->c[s - 1] != 1.) scheme->stiffly_accurate = PETSC_FALSE; 273 for (j = 0; j < s; j++) 274 if (a[(s - 1) * s + j] != b[j]) scheme->stiffly_accurate = PETSC_FALSE; 275 for (j = 0; j < r; j++) 276 if (u[(s - 1) * r + j] != v[j]) scheme->stiffly_accurate = PETSC_FALSE; 277 scheme->fsal = scheme->stiffly_accurate; /* FSAL is stronger */ 278 for (j = 0; j < s - 1; j++) 279 if (r > 1 && b[1 * s + j] != 0.) scheme->fsal = PETSC_FALSE; 280 if (b[1 * s + r - 1] != 1.) scheme->fsal = PETSC_FALSE; 281 for (j = 0; j < r; j++) 282 if (r > 1 && v[1 * r + j] != 0.) scheme->fsal = PETSC_FALSE; 283 284 *inscheme = scheme; 285 PetscFunctionReturn(PETSC_SUCCESS); 286 } 287 288 static PetscErrorCode TSGLLESchemeDestroy(TSGLLEScheme sc) 289 { 290 PetscFunctionBegin; 291 PetscCall(PetscFree5(sc->c, sc->a, sc->b, sc->u, sc->v)); 292 PetscCall(PetscFree6(sc->alpha, sc->beta, sc->gamma, sc->phi, sc->psi, sc->stage_error)); 293 PetscCall(PetscFree(sc)); 294 PetscFunctionReturn(PETSC_SUCCESS); 295 } 296 297 static PetscErrorCode TSGLLEDestroy_Default(TS_GLLE *gl) 298 { 299 PetscInt i; 300 301 PetscFunctionBegin; 302 for (i = 0; i < gl->nschemes; i++) { 303 if (gl->schemes[i]) PetscCall(TSGLLESchemeDestroy(gl->schemes[i])); 304 } 305 PetscCall(PetscFree(gl->schemes)); 306 gl->nschemes = 0; 307 PetscCall(PetscMemzero(gl->type_name, sizeof(gl->type_name))); 308 PetscFunctionReturn(PETSC_SUCCESS); 309 } 310 311 static PetscErrorCode TSGLLEViewTable_Private(PetscViewer viewer, PetscInt m, PetscInt n, const PetscScalar a[], const char name[]) 312 { 313 PetscBool iascii; 314 PetscInt i, j; 315 316 PetscFunctionBegin; 317 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii)); 318 if (iascii) { 319 PetscCall(PetscViewerASCIIPrintf(viewer, "%30s = [", name)); 320 for (i = 0; i < m; i++) { 321 if (i) PetscCall(PetscViewerASCIIPrintf(viewer, "%30s [", "")); 322 PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_FALSE)); 323 for (j = 0; j < n; j++) PetscCall(PetscViewerASCIIPrintf(viewer, " %12.8g", (double)PetscRealPart(a[i * n + j]))); 324 PetscCall(PetscViewerASCIIPrintf(viewer, "]\n")); 325 PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_TRUE)); 326 } 327 } 328 PetscFunctionReturn(PETSC_SUCCESS); 329 } 330 331 static PetscErrorCode TSGLLESchemeView(TSGLLEScheme sc, PetscBool view_details, PetscViewer viewer) 332 { 333 PetscBool iascii; 334 335 PetscFunctionBegin; 336 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii)); 337 if (iascii) { 338 PetscCall(PetscViewerASCIIPrintf(viewer, "GL scheme p,q,r,s = %" PetscInt_FMT ",%" PetscInt_FMT ",%" PetscInt_FMT ",%" PetscInt_FMT "\n", sc->p, sc->q, sc->r, sc->s)); 339 PetscCall(PetscViewerASCIIPushTab(viewer)); 340 PetscCall(PetscViewerASCIIPrintf(viewer, "Stiffly accurate: %s, FSAL: %s\n", sc->stiffly_accurate ? "yes" : "no", sc->fsal ? "yes" : "no")); 341 PetscCall(PetscViewerASCIIPrintf(viewer, "Leading error constants: %10.3e %10.3e %10.3e\n", (double)PetscRealPart(sc->alpha[0]), (double)PetscRealPart(sc->beta[0]), (double)PetscRealPart(sc->gamma[0]))); 342 PetscCall(TSGLLEViewTable_Private(viewer, 1, sc->s, sc->c, "Abscissas c")); 343 if (view_details) { 344 PetscCall(TSGLLEViewTable_Private(viewer, sc->s, sc->s, sc->a, "A")); 345 PetscCall(TSGLLEViewTable_Private(viewer, sc->r, sc->s, sc->b, "B")); 346 PetscCall(TSGLLEViewTable_Private(viewer, sc->s, sc->r, sc->u, "U")); 347 PetscCall(TSGLLEViewTable_Private(viewer, sc->r, sc->r, sc->v, "V")); 348 349 PetscCall(TSGLLEViewTable_Private(viewer, 3, sc->s, sc->phi, "Error estimate phi")); 350 PetscCall(TSGLLEViewTable_Private(viewer, 3, sc->r, sc->psi, "Error estimate psi")); 351 PetscCall(TSGLLEViewTable_Private(viewer, 1, sc->r, sc->alpha, "Modify alpha")); 352 PetscCall(TSGLLEViewTable_Private(viewer, 1, sc->r, sc->beta, "Modify beta")); 353 PetscCall(TSGLLEViewTable_Private(viewer, 1, sc->r, sc->gamma, "Modify gamma")); 354 PetscCall(TSGLLEViewTable_Private(viewer, 1, sc->s, sc->stage_error, "Stage error xi")); 355 } 356 PetscCall(PetscViewerASCIIPopTab(viewer)); 357 } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Viewer type %s not supported", ((PetscObject)viewer)->type_name); 358 PetscFunctionReturn(PETSC_SUCCESS); 359 } 360 361 static PetscErrorCode TSGLLEEstimateHigherMoments_Default(TSGLLEScheme sc, PetscReal h, Vec Ydot[], Vec Xold[], Vec hm[]) 362 { 363 PetscInt i; 364 365 PetscFunctionBegin; 366 PetscCheck(sc->r <= 64 && sc->s <= 64, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Ridiculous number of stages or items passed between stages"); 367 /* build error vectors*/ 368 for (i = 0; i < 3; i++) { 369 PetscScalar phih[64]; 370 PetscInt j; 371 for (j = 0; j < sc->s; j++) phih[j] = sc->phi[i * sc->s + j] * h; 372 PetscCall(VecZeroEntries(hm[i])); 373 PetscCall(VecMAXPY(hm[i], sc->s, phih, Ydot)); 374 PetscCall(VecMAXPY(hm[i], sc->r, &sc->psi[i * sc->r], Xold)); 375 } 376 PetscFunctionReturn(PETSC_SUCCESS); 377 } 378 379 static PetscErrorCode TSGLLECompleteStep_Rescale(TSGLLEScheme sc, PetscReal h, TSGLLEScheme next_sc, PetscReal next_h, Vec Ydot[], Vec Xold[], Vec X[]) 380 { 381 PetscScalar brow[32], vrow[32]; 382 PetscInt i, j, r, s; 383 384 PetscFunctionBegin; 385 /* Build the new solution from (X,Ydot) */ 386 r = sc->r; 387 s = sc->s; 388 for (i = 0; i < r; i++) { 389 PetscCall(VecZeroEntries(X[i])); 390 for (j = 0; j < s; j++) brow[j] = h * sc->b[i * s + j]; 391 PetscCall(VecMAXPY(X[i], s, brow, Ydot)); 392 for (j = 0; j < r; j++) vrow[j] = sc->v[i * r + j]; 393 PetscCall(VecMAXPY(X[i], r, vrow, Xold)); 394 } 395 PetscFunctionReturn(PETSC_SUCCESS); 396 } 397 398 static PetscErrorCode TSGLLECompleteStep_RescaleAndModify(TSGLLEScheme sc, PetscReal h, TSGLLEScheme next_sc, PetscReal next_h, Vec Ydot[], Vec Xold[], Vec X[]) 399 { 400 PetscScalar brow[32], vrow[32]; 401 PetscReal ratio; 402 PetscInt i, j, p, r, s; 403 404 PetscFunctionBegin; 405 /* Build the new solution from (X,Ydot) */ 406 p = sc->p; 407 r = sc->r; 408 s = sc->s; 409 ratio = next_h / h; 410 for (i = 0; i < r; i++) { 411 PetscCall(VecZeroEntries(X[i])); 412 for (j = 0; j < s; j++) { 413 brow[j] = h * (PetscPowRealInt(ratio, i) * sc->b[i * s + j] + (PetscPowRealInt(ratio, i) - PetscPowRealInt(ratio, p + 1)) * (+sc->alpha[i] * sc->phi[0 * s + j]) + (PetscPowRealInt(ratio, i) - PetscPowRealInt(ratio, p + 2)) * (+sc->beta[i] * sc->phi[1 * s + j] + sc->gamma[i] * sc->phi[2 * s + j])); 414 } 415 PetscCall(VecMAXPY(X[i], s, brow, Ydot)); 416 for (j = 0; j < r; j++) { 417 vrow[j] = (PetscPowRealInt(ratio, i) * sc->v[i * r + j] + (PetscPowRealInt(ratio, i) - PetscPowRealInt(ratio, p + 1)) * (+sc->alpha[i] * sc->psi[0 * r + j]) + (PetscPowRealInt(ratio, i) - PetscPowRealInt(ratio, p + 2)) * (+sc->beta[i] * sc->psi[1 * r + j] + sc->gamma[i] * sc->psi[2 * r + j])); 418 } 419 PetscCall(VecMAXPY(X[i], r, vrow, Xold)); 420 } 421 if (r < next_sc->r) { 422 PetscCheck(r + 1 == next_sc->r, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Cannot accommodate jump in r greater than 1"); 423 PetscCall(VecZeroEntries(X[r])); 424 for (j = 0; j < s; j++) brow[j] = h * PetscPowRealInt(ratio, p + 1) * sc->phi[0 * s + j]; 425 PetscCall(VecMAXPY(X[r], s, brow, Ydot)); 426 for (j = 0; j < r; j++) vrow[j] = PetscPowRealInt(ratio, p + 1) * sc->psi[0 * r + j]; 427 PetscCall(VecMAXPY(X[r], r, vrow, Xold)); 428 } 429 PetscFunctionReturn(PETSC_SUCCESS); 430 } 431 432 static PetscErrorCode TSGLLECreate_IRKS(TS ts) 433 { 434 TS_GLLE *gl = (TS_GLLE *)ts->data; 435 436 PetscFunctionBegin; 437 gl->Destroy = TSGLLEDestroy_Default; 438 gl->EstimateHigherMoments = TSGLLEEstimateHigherMoments_Default; 439 gl->CompleteStep = TSGLLECompleteStep_RescaleAndModify; 440 PetscCall(PetscMalloc1(10, &gl->schemes)); 441 gl->nschemes = 0; 442 443 { 444 /* p=1,q=1, r=s=2, A- and L-stable with error estimates of order 2 and 3 445 * Listed in Butcher & Podhaisky 2006. On error estimation in general linear methods for stiff ODE. 446 * irks(0.3,0,[.3,1],[1],1) 447 * Note: can be made to have classical order (not stage order) 2 by replacing 0.3 with 1-sqrt(1/2) 448 * but doing so would sacrifice the error estimator. 449 */ 450 const PetscScalar c[2] = {3. / 10., 1.}; 451 const PetscScalar a[2][2] = { 452 {3. / 10., 0 }, 453 {7. / 10., 3. / 10.} 454 }; 455 const PetscScalar b[2][2] = { 456 {7. / 10., 3. / 10.}, 457 {0, 1 } 458 }; 459 const PetscScalar u[2][2] = { 460 {1, 0}, 461 {1, 0} 462 }; 463 const PetscScalar v[2][2] = { 464 {1, 0}, 465 {0, 0} 466 }; 467 PetscCall(TSGLLESchemeCreate(1, 1, 2, 2, c, *a, *b, *u, *v, &gl->schemes[gl->nschemes++])); 468 } 469 470 { 471 /* p=q=2, r=s=3: irks(4/9,0,[1:3]/3,[0.33852],1) */ 472 /* http://www.math.auckland.ac.nz/~hpod/atlas/i2a.html */ 473 const PetscScalar c[3] = {1. / 3., 2. / 3., 1}; 474 const PetscScalar a[3][3] = { 475 {4. / 9., 0, 0 }, 476 {1.03750643704090e+00, 4. / 9., 0 }, 477 {7.67024779410304e-01, -3.81140216918943e-01, 4. / 9.} 478 }; 479 const PetscScalar b[3][3] = { 480 {0.767024779410304, -0.381140216918943, 4. / 9. }, 481 {0.000000000000000, 0.000000000000000, 1.000000000000000}, 482 {-2.075048385225385, 0.621728385225383, 1.277197204924873} 483 }; 484 const PetscScalar u[3][3] = { 485 {1.0000000000000000, -0.1111111111111109, -0.0925925925925922}, 486 {1.0000000000000000, -0.8152842148186744, -0.4199095530877056}, 487 {1.0000000000000000, 0.1696709930641948, 0.0539741070314165 } 488 }; 489 const PetscScalar v[3][3] = { 490 {1.0000000000000000, 0.1696709930641948, 0.0539741070314165}, 491 {0.000000000000000, 0.000000000000000, 0.000000000000000 }, 492 {0.000000000000000, 0.176122795075129, 0.000000000000000 } 493 }; 494 PetscCall(TSGLLESchemeCreate(2, 2, 3, 3, c, *a, *b, *u, *v, &gl->schemes[gl->nschemes++])); 495 } 496 { 497 /* p=q=3, r=s=4: irks(9/40,0,[1:4]/4,[0.3312 1.0050],[0.49541 1;1 0]) */ 498 const PetscScalar c[4] = {0.25, 0.5, 0.75, 1.0}; 499 const PetscScalar a[4][4] = { 500 {9. / 40., 0, 0, 0 }, 501 {2.11286958887701e-01, 9. / 40., 0, 0 }, 502 {9.46338294287584e-01, -3.42942861246094e-01, 9. / 40., 0 }, 503 {0.521490453970721, -0.662474225622980, 0.490476425459734, 9. / 40.} 504 }; 505 const PetscScalar b[4][4] = { 506 {0.521490453970721, -0.662474225622980, 0.490476425459734, 9. / 40. }, 507 {0.000000000000000, 0.000000000000000, 0.000000000000000, 1.000000000000000}, 508 {-0.084677029310348, 1.390757514776085, -1.568157386206001, 2.023192696767826}, 509 {0.465383797936408, 1.478273530625148, -1.930836081010182, 1.644872111193354} 510 }; 511 const PetscScalar u[4][4] = { 512 {1.00000000000000000, 0.02500000000001035, -0.02499999999999053, -0.00442708333332865}, 513 {1.00000000000000000, 0.06371304111232945, -0.04032173972189845, -0.01389438413189452}, 514 {1.00000000000000000, -0.07839543304147778, 0.04738685705116663, 0.02032603595928376 }, 515 {1.00000000000000000, 0.42550734619251651, 0.10800718022400080, -0.01726712647760034} 516 }; 517 const PetscScalar v[4][4] = { 518 {1.00000000000000000, 0.42550734619251651, 0.10800718022400080, -0.01726712647760034}, 519 {0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000 }, 520 {0.000000000000000, -1.761115796027561, -0.521284157173780, 0.258249384305463 }, 521 {0.000000000000000, -1.657693358744728, -1.052227765232394, 0.521284157173780 } 522 }; 523 PetscCall(TSGLLESchemeCreate(3, 3, 4, 4, c, *a, *b, *u, *v, &gl->schemes[gl->nschemes++])); 524 } 525 { 526 /* p=q=4, r=s=5: 527 irks(3/11,0,[1:5]/5, [0.1715 -0.1238 0.6617],... 528 [ -0.0812 0.4079 1.0000 529 1.0000 0 0 530 0.8270 1.0000 0]) 531 */ 532 const PetscScalar c[5] = {0.2, 0.4, 0.6, 0.8, 1.0}; 533 const PetscScalar a[5][5] = { 534 {2.72727272727352e-01, 0.00000000000000e+00, 0.00000000000000e+00, 0.00000000000000e+00, 0.00000000000000e+00}, 535 {-1.03980153733431e-01, 2.72727272727405e-01, 0.00000000000000e+00, 0.00000000000000e+00, 0.00000000000000e+00}, 536 {-1.58615400341492e+00, 7.44168951881122e-01, 2.72727272727309e-01, 0.00000000000000e+00, 0.00000000000000e+00}, 537 {-8.73658042865628e-01, 5.37884671894595e-01, -1.63298538799523e-01, 2.72727272726996e-01, 0.00000000000000e+00}, 538 {2.95489397443992e-01, -1.18481693910097e+00, -6.68029812659953e-01, 1.00716687860943e+00, 2.72727272727288e-01} 539 }; 540 const PetscScalar b[5][5] = { 541 {2.95489397443992e-01, -1.18481693910097e+00, -6.68029812659953e-01, 1.00716687860943e+00, 2.72727272727288e-01}, 542 {0.00000000000000e+00, 1.11022302462516e-16, -2.22044604925031e-16, 0.00000000000000e+00, 1.00000000000000e+00}, 543 {-4.05882503986005e+00, -4.00924006567769e+00, -1.38930610972481e+00, 4.45223930308488e+00, 6.32331093108427e-01}, 544 {8.35690179937017e+00, -2.26640927349732e+00, 6.86647884973826e+00, -5.22595158025740e+00, 4.50893068837431e+00}, 545 {1.27656267027479e+01, 2.80882153840821e+00, 8.91173096522890e+00, -1.07936444078906e+01, 4.82534148988854e+00} 546 }; 547 const PetscScalar u[5][5] = { 548 {1.00000000000000e+00, -7.27272727273551e-02, -3.45454545454419e-02, -4.12121212119565e-03, -2.96969696964014e-04}, 549 {1.00000000000000e+00, 2.31252881006154e-01, -8.29487834416481e-03, -9.07191207681020e-03, -1.70378403743473e-03}, 550 {1.00000000000000e+00, 1.16925777880663e+00, 3.59268562942635e-02, -4.09013451730615e-02, -1.02411119670164e-02}, 551 {1.00000000000000e+00, 1.02634463704356e+00, 1.59375044913405e-01, 1.89673015035370e-03, -4.89987231897569e-03}, 552 {1.00000000000000e+00, 1.27746320298021e+00, 2.37186008132728e-01, -8.28694373940065e-02, -5.34396510196430e-02} 553 }; 554 const PetscScalar v[5][5] = { 555 {1.00000000000000e+00, 1.27746320298021e+00, 2.37186008132728e-01, -8.28694373940065e-02, -5.34396510196430e-02}, 556 {0.00000000000000e+00, -1.77635683940025e-15, -1.99840144432528e-15, -9.99200722162641e-16, -3.33066907387547e-16}, 557 {0.00000000000000e+00, 4.37280081906924e+00, 5.49221645016377e-02, -8.88913177394943e-02, 1.12879077989154e-01 }, 558 {0.00000000000000e+00, -1.22399504837280e+01, -5.21287338448645e+00, -8.03952325565291e-01, 4.60298678047147e-01 }, 559 {0.00000000000000e+00, -1.85178762883829e+01, -5.21411849862624e+00, -1.04283436528809e+00, 7.49030161063651e-01 } 560 }; 561 PetscCall(TSGLLESchemeCreate(4, 4, 5, 5, c, *a, *b, *u, *v, &gl->schemes[gl->nschemes++])); 562 } 563 { 564 /* p=q=5, r=s=6; 565 irks(1/3,0,[1:6]/6,... 566 [-0.0489 0.4228 -0.8814 0.9021],... 567 [-0.3474 -0.6617 0.6294 0.2129 568 0.0044 -0.4256 -0.1427 -0.8936 569 -0.8267 0.4821 0.1371 -0.2557 570 -0.4426 -0.3855 -0.7514 0.3014]) 571 */ 572 const PetscScalar c[6] = {1. / 6, 2. / 6, 3. / 6, 4. / 6, 5. / 6, 1.}; 573 const PetscScalar a[6][6] = { 574 {3.33333333333940e-01, 0, 0, 0, 0, 0 }, 575 {-8.64423857333350e-02, 3.33333333332888e-01, 0, 0, 0, 0 }, 576 {-2.16850174258252e+00, -2.23619072028839e+00, 3.33333333335204e-01, 0, 0, 0 }, 577 {-4.73160970138997e+00, -3.89265344629268e+00, -2.76318716520933e-01, 3.33333333335759e-01, 0, 0 }, 578 {-6.75187540297338e+00, -7.90756533769377e+00, 7.90245051802259e-01, -4.48352364517632e-01, 3.33333333328483e-01, 0 }, 579 {-4.26488287921548e+00, -1.19320395589302e+01, 3.38924509887755e+00, -2.23969848002481e+00, 6.62807710124007e-01, 3.33333333335440e-01} 580 }; 581 const PetscScalar b[6][6] = { 582 {-4.26488287921548e+00, -1.19320395589302e+01, 3.38924509887755e+00, -2.23969848002481e+00, 6.62807710124007e-01, 3.33333333335440e-01 }, 583 {-8.88178419700125e-16, 4.44089209850063e-16, -1.54737334057131e-15, -8.88178419700125e-16, 0.00000000000000e+00, 1.00000000000001e+00 }, 584 {-2.87780425770651e+01, -1.13520448264971e+01, 2.62002318943161e+01, 2.56943874812797e+01, -3.06702268304488e+01, 6.68067773510103e+00 }, 585 {5.47971245256474e+01, 6.80366875868284e+01, -6.50952588861999e+01, -8.28643975339097e+01, 8.17416943896414e+01, -1.17819043489036e+01}, 586 {-2.33332114788869e+02, 6.12942539462634e+01, -4.91850135865944e+01, 1.82716844135480e+02, -1.29788173979395e+02, 3.09968095651099e+01 }, 587 {-1.72049132343751e+02, 8.60194713593999e+00, 7.98154219170200e-01, 1.50371386053218e+02, -1.18515423962066e+02, 2.50898277784663e+01 } 588 }; 589 const PetscScalar u[6][6] = { 590 {1.00000000000000e+00, -1.66666666666870e-01, -4.16666666664335e-02, -3.85802469124815e-03, -2.25051440302250e-04, -9.64506172339142e-06}, 591 {1.00000000000000e+00, 8.64423857327162e-02, -4.11484912671353e-02, -1.11450903217645e-02, -1.47651050487126e-03, -1.34395070766826e-04}, 592 {1.00000000000000e+00, 4.57135912953434e+00, 1.06514719719137e+00, 1.33517564218007e-01, 1.11365952968659e-02, 6.12382756769504e-04 }, 593 {1.00000000000000e+00, 9.23391519753404e+00, 2.22431212392095e+00, 2.91823807741891e-01, 2.52058456411084e-02, 1.22800542949647e-03 }, 594 {1.00000000000000e+00, 1.48175480533865e+01, 3.73439117461835e+00, 5.14648336541804e-01, 4.76430038853402e-02, 2.56798515502156e-03 }, 595 {1.00000000000000e+00, 1.50512347758335e+01, 4.10099701165164e+00, 5.66039141003603e-01, 3.91213893800891e-02, -2.99136269067853e-03} 596 }; 597 const PetscScalar v[6][6] = { 598 {1.00000000000000e+00, 1.50512347758335e+01, 4.10099701165164e+00, 5.66039141003603e-01, 3.91213893800891e-02, -2.99136269067853e-03}, 599 {0.00000000000000e+00, -4.88498130835069e-15, -6.43929354282591e-15, -3.55271367880050e-15, -1.22124532708767e-15, -3.12250225675825e-16}, 600 {0.00000000000000e+00, 1.22250171233141e+01, -1.77150760606169e+00, 3.54516769879390e-01, 6.22298845883398e-01, 2.31647447450276e-01 }, 601 {0.00000000000000e+00, -4.48339457331040e+01, -3.57363126641880e-01, 5.18750173123425e-01, 6.55727990241799e-02, 1.63175368287079e-01 }, 602 {0.00000000000000e+00, 1.37297394708005e+02, -1.60145272991317e+00, -5.05319555199441e+00, 1.55328940390990e-01, 9.16629423682464e-01 }, 603 {0.00000000000000e+00, 1.05703241119022e+02, -1.16610260983038e+00, -2.99767252773859e+00, -1.13472315553890e-01, 1.09742849254729e+00 } 604 }; 605 PetscCall(TSGLLESchemeCreate(5, 5, 6, 6, c, *a, *b, *u, *v, &gl->schemes[gl->nschemes++])); 606 } 607 PetscFunctionReturn(PETSC_SUCCESS); 608 } 609 610 /*@ 611 TSGLLESetType - sets the class of general linear method, `TSGLLE` to use for time-stepping 612 613 Collective 614 615 Input Parameters: 616 + ts - the `TS` context 617 - type - a method 618 619 Options Database Key: 620 . -ts_gl_type <type> - sets the method, use -help for a list of available method (e.g. irks) 621 622 Level: intermediate 623 624 Notes: 625 See "petsc/include/petscts.h" for available methods (for instance) 626 . TSGLLE_IRKS - Diagonally implicit methods with inherent Runge-Kutta stability (for stiff problems) 627 628 Normally, it is best to use the `TSSetFromOptions()` command and then set the `TSGLLE` type 629 from the options database rather than by using this routine. Using the options database 630 provides the user with maximum flexibility in evaluating the many different solvers. The 631 `TSGLLESetType()` routine is provided for those situations where it is necessary to set the 632 timestepping solver independently of the command line or options database. This might be the 633 case, for example, when the choice of solver changes during the execution of the program, and 634 the user's application is taking responsibility for choosing the appropriate method. 635 636 .seealso: [](ch_ts), `TS`, `TSGLLEType`, `TSGLLE` 637 @*/ 638 PetscErrorCode TSGLLESetType(TS ts, TSGLLEType type) 639 { 640 PetscFunctionBegin; 641 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 642 PetscAssertPointer(type, 2); 643 PetscTryMethod(ts, "TSGLLESetType_C", (TS, TSGLLEType), (ts, type)); 644 PetscFunctionReturn(PETSC_SUCCESS); 645 } 646 647 /*@C 648 TSGLLESetAcceptType - sets the acceptance test for `TSGLLE` 649 650 Logically Collective 651 652 Input Parameters: 653 + ts - the `TS` context 654 - type - the type 655 656 Options Database Key: 657 . -ts_gl_accept_type <type> - sets the method used to determine whether to accept or reject a step 658 659 Level: intermediate 660 661 Notes: 662 Time integrators that need to control error must have the option to reject a time step based 663 on local error estimates. This function allows different schemes to be set. 664 665 .seealso: [](ch_ts), `TS`, `TSGLLE`, `TSGLLEAcceptRegister()`, `TSGLLEAdapt` 666 @*/ 667 PetscErrorCode TSGLLESetAcceptType(TS ts, TSGLLEAcceptType type) 668 { 669 PetscFunctionBegin; 670 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 671 PetscAssertPointer(type, 2); 672 PetscTryMethod(ts, "TSGLLESetAcceptType_C", (TS, TSGLLEAcceptType), (ts, type)); 673 PetscFunctionReturn(PETSC_SUCCESS); 674 } 675 676 /*@ 677 TSGLLEGetAdapt - gets the `TSGLLEAdapt` object from the `TS` 678 679 Not Collective 680 681 Input Parameter: 682 . ts - the `TS` context 683 684 Output Parameter: 685 . adapt - the `TSGLLEAdapt` context 686 687 Level: advanced 688 689 Note: 690 This allows the user set options on the `TSGLLEAdapt` object. Usually it is better to do this 691 using the options database, so this function is rarely needed. 692 693 .seealso: [](ch_ts), `TS`, `TSGLLE`, `TSGLLEAdapt`, `TSGLLEAdaptRegister()` 694 @*/ 695 PetscErrorCode TSGLLEGetAdapt(TS ts, TSGLLEAdapt *adapt) 696 { 697 PetscFunctionBegin; 698 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 699 PetscAssertPointer(adapt, 2); 700 PetscUseMethod(ts, "TSGLLEGetAdapt_C", (TS, TSGLLEAdapt *), (ts, adapt)); 701 PetscFunctionReturn(PETSC_SUCCESS); 702 } 703 704 static PetscErrorCode TSGLLEAccept_Always(TS ts, PetscReal tleft, PetscReal h, const PetscReal enorms[], PetscBool *accept) 705 { 706 PetscFunctionBegin; 707 *accept = PETSC_TRUE; 708 PetscFunctionReturn(PETSC_SUCCESS); 709 } 710 711 static PetscErrorCode TSGLLEUpdateWRMS(TS ts) 712 { 713 TS_GLLE *gl = (TS_GLLE *)ts->data; 714 PetscScalar *x, *w; 715 PetscInt n, i; 716 717 PetscFunctionBegin; 718 PetscCall(VecGetArray(gl->X[0], &x)); 719 PetscCall(VecGetArray(gl->W, &w)); 720 PetscCall(VecGetLocalSize(gl->W, &n)); 721 for (i = 0; i < n; i++) w[i] = 1. / (gl->wrms_atol + gl->wrms_rtol * PetscAbsScalar(x[i])); 722 PetscCall(VecRestoreArray(gl->X[0], &x)); 723 PetscCall(VecRestoreArray(gl->W, &w)); 724 PetscFunctionReturn(PETSC_SUCCESS); 725 } 726 727 static PetscErrorCode TSGLLEVecNormWRMS(TS ts, Vec X, PetscReal *nrm) 728 { 729 TS_GLLE *gl = (TS_GLLE *)ts->data; 730 PetscScalar *x, *w; 731 PetscReal sum = 0.0, gsum; 732 PetscInt n, N, i; 733 734 PetscFunctionBegin; 735 PetscCall(VecGetArray(X, &x)); 736 PetscCall(VecGetArray(gl->W, &w)); 737 PetscCall(VecGetLocalSize(gl->W, &n)); 738 for (i = 0; i < n; i++) sum += PetscAbsScalar(PetscSqr(x[i] * w[i])); 739 PetscCall(VecRestoreArray(X, &x)); 740 PetscCall(VecRestoreArray(gl->W, &w)); 741 PetscCallMPI(MPIU_Allreduce(&sum, &gsum, 1, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)ts))); 742 PetscCall(VecGetSize(gl->W, &N)); 743 *nrm = PetscSqrtReal(gsum / (1. * N)); 744 PetscFunctionReturn(PETSC_SUCCESS); 745 } 746 747 static PetscErrorCode TSGLLESetType_GLLE(TS ts, TSGLLEType type) 748 { 749 PetscBool same; 750 TS_GLLE *gl = (TS_GLLE *)ts->data; 751 PetscErrorCode (*r)(TS); 752 753 PetscFunctionBegin; 754 if (gl->type_name[0]) { 755 PetscCall(PetscStrcmp(gl->type_name, type, &same)); 756 if (same) PetscFunctionReturn(PETSC_SUCCESS); 757 PetscCall((*gl->Destroy)(gl)); 758 } 759 760 PetscCall(PetscFunctionListFind(TSGLLEList, type, &r)); 761 PetscCheck(r, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown TSGLLE type \"%s\" given", type); 762 PetscCall((*r)(ts)); 763 PetscCall(PetscStrncpy(gl->type_name, type, sizeof(gl->type_name))); 764 PetscFunctionReturn(PETSC_SUCCESS); 765 } 766 767 static PetscErrorCode TSGLLESetAcceptType_GLLE(TS ts, TSGLLEAcceptType type) 768 { 769 TSGLLEAcceptFn *r; 770 TS_GLLE *gl = (TS_GLLE *)ts->data; 771 772 PetscFunctionBegin; 773 PetscCall(PetscFunctionListFind(TSGLLEAcceptList, type, &r)); 774 PetscCheck(r, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown TSGLLEAccept type \"%s\" given", type); 775 gl->Accept = r; 776 PetscCall(PetscStrncpy(gl->accept_name, type, sizeof(gl->accept_name))); 777 PetscFunctionReturn(PETSC_SUCCESS); 778 } 779 780 static PetscErrorCode TSGLLEGetAdapt_GLLE(TS ts, TSGLLEAdapt *adapt) 781 { 782 TS_GLLE *gl = (TS_GLLE *)ts->data; 783 784 PetscFunctionBegin; 785 if (!gl->adapt) { 786 PetscCall(TSGLLEAdaptCreate(PetscObjectComm((PetscObject)ts), &gl->adapt)); 787 PetscCall(PetscObjectIncrementTabLevel((PetscObject)gl->adapt, (PetscObject)ts, 1)); 788 } 789 *adapt = gl->adapt; 790 PetscFunctionReturn(PETSC_SUCCESS); 791 } 792 793 static PetscErrorCode TSGLLEChooseNextScheme(TS ts, PetscReal h, const PetscReal hmnorm[], PetscInt *next_scheme, PetscReal *next_h, PetscBool *finish) 794 { 795 TS_GLLE *gl = (TS_GLLE *)ts->data; 796 PetscInt i, n, cur_p, cur, next_sc, candidates[64], orders[64]; 797 PetscReal errors[64], costs[64], tleft; 798 799 PetscFunctionBegin; 800 cur = -1; 801 cur_p = gl->schemes[gl->current_scheme]->p; 802 tleft = ts->max_time - (ts->ptime + ts->time_step); 803 for (i = 0, n = 0; i < gl->nschemes; i++) { 804 TSGLLEScheme sc = gl->schemes[i]; 805 if (sc->p < gl->min_order || gl->max_order < sc->p) continue; 806 if (sc->p == cur_p - 1) errors[n] = PetscAbsScalar(sc->alpha[0]) * hmnorm[0]; 807 else if (sc->p == cur_p) errors[n] = PetscAbsScalar(sc->alpha[0]) * hmnorm[1]; 808 else if (sc->p == cur_p + 1) errors[n] = PetscAbsScalar(sc->alpha[0]) * (hmnorm[2] + hmnorm[3]); 809 else continue; 810 candidates[n] = i; 811 orders[n] = PetscMin(sc->p, sc->q); /* order of global truncation error */ 812 costs[n] = sc->s; /* estimate the cost as the number of stages */ 813 if (i == gl->current_scheme) cur = n; 814 n++; 815 } 816 PetscCheck(cur >= 0 && gl->nschemes > cur, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Current scheme not found in scheme list"); 817 PetscCall(TSGLLEAdaptChoose(gl->adapt, n, orders, errors, costs, cur, h, tleft, &next_sc, next_h, finish)); 818 *next_scheme = candidates[next_sc]; 819 PetscCall(PetscInfo(ts, "Adapt chose scheme %" PetscInt_FMT " (%" PetscInt_FMT ",%" PetscInt_FMT ",%" PetscInt_FMT ",%" PetscInt_FMT ") with step size %6.2e, finish=%s\n", *next_scheme, gl->schemes[*next_scheme]->p, gl->schemes[*next_scheme]->q, 820 gl->schemes[*next_scheme]->r, gl->schemes[*next_scheme]->s, (double)*next_h, PetscBools[*finish])); 821 PetscFunctionReturn(PETSC_SUCCESS); 822 } 823 824 static PetscErrorCode TSGLLEGetMaxSizes(TS ts, PetscInt *max_r, PetscInt *max_s) 825 { 826 TS_GLLE *gl = (TS_GLLE *)ts->data; 827 828 PetscFunctionBegin; 829 *max_r = gl->schemes[gl->nschemes - 1]->r; 830 *max_s = gl->schemes[gl->nschemes - 1]->s; 831 PetscFunctionReturn(PETSC_SUCCESS); 832 } 833 834 static PetscErrorCode TSSolve_GLLE(TS ts) 835 { 836 TS_GLLE *gl = (TS_GLLE *)ts->data; 837 PetscInt i, k, its, lits, max_r, max_s; 838 PetscBool final_step, finish; 839 SNESConvergedReason snesreason; 840 841 PetscFunctionBegin; 842 PetscCall(TSMonitor(ts, ts->steps, ts->ptime, ts->vec_sol)); 843 844 PetscCall(TSGLLEGetMaxSizes(ts, &max_r, &max_s)); 845 PetscCall(VecCopy(ts->vec_sol, gl->X[0])); 846 for (i = 1; i < max_r; i++) PetscCall(VecZeroEntries(gl->X[i])); 847 PetscCall(TSGLLEUpdateWRMS(ts)); 848 849 if (0) { 850 /* Find consistent initial data for DAE */ 851 gl->stage_time = ts->ptime + ts->time_step; 852 gl->scoeff = 1.; 853 gl->stage = 0; 854 855 PetscCall(VecCopy(ts->vec_sol, gl->Z)); 856 PetscCall(VecCopy(ts->vec_sol, gl->Y)); 857 PetscCall(SNESSolve(ts->snes, NULL, gl->Y)); 858 PetscCall(SNESGetIterationNumber(ts->snes, &its)); 859 PetscCall(SNESGetLinearSolveIterations(ts->snes, &lits)); 860 PetscCall(SNESGetConvergedReason(ts->snes, &snesreason)); 861 862 ts->snes_its += its; 863 ts->ksp_its += lits; 864 if (snesreason < 0 && ts->max_snes_failures != PETSC_UNLIMITED && ++ts->num_snes_failures >= ts->max_snes_failures) { 865 ts->reason = TS_DIVERGED_NONLINEAR_SOLVE; 866 PetscCall(PetscInfo(ts, "Step=%" PetscInt_FMT ", nonlinear solve failures %" PetscInt_FMT " greater than current TS allowed, stopping solve\n", ts->steps, ts->num_snes_failures)); 867 PetscFunctionReturn(PETSC_SUCCESS); 868 } 869 } 870 871 PetscCheck(gl->current_scheme >= 0, PETSC_COMM_SELF, PETSC_ERR_ORDER, "A starting scheme has not been provided"); 872 873 for (k = 0, final_step = PETSC_FALSE, finish = PETSC_FALSE; k < ts->max_steps && !finish; k++) { 874 PetscInt j, r, s, next_scheme = 0, rejections; 875 PetscReal h, hmnorm[4], enorm[3], next_h; 876 PetscBool accept; 877 const PetscScalar *c, *a, *u; 878 Vec *X, *Ydot, Y; 879 TSGLLEScheme scheme = gl->schemes[gl->current_scheme]; 880 881 r = scheme->r; 882 s = scheme->s; 883 c = scheme->c; 884 a = scheme->a; 885 u = scheme->u; 886 h = ts->time_step; 887 X = gl->X; 888 Ydot = gl->Ydot; 889 Y = gl->Y; 890 891 if (ts->ptime > ts->max_time) break; 892 893 /* 894 We only call PreStep at the start of each STEP, not each STAGE. This is because it is 895 possible to fail (have to restart a step) after multiple stages. 896 */ 897 PetscCall(TSPreStep(ts)); 898 899 rejections = 0; 900 while (1) { 901 for (i = 0; i < s; i++) { 902 PetscScalar shift; 903 gl->scoeff = 1. / PetscRealPart(a[i * s + i]); 904 shift = gl->scoeff / ts->time_step; 905 gl->stage = i; 906 gl->stage_time = ts->ptime + PetscRealPart(c[i]) * h; 907 908 /* 909 * Stage equation: Y = h A Y' + U X 910 * We assume that A is lower-triangular so that we can solve the stages (Y,Y') sequentially 911 * Build the affine vector z_i = -[1/(h a_ii)](h sum_j a_ij y'_j + sum_j u_ij x_j) 912 * Then y'_i = z + 1/(h a_ii) y_i 913 */ 914 PetscCall(VecZeroEntries(gl->Z)); 915 for (j = 0; j < r; j++) PetscCall(VecAXPY(gl->Z, -shift * u[i * r + j], X[j])); 916 for (j = 0; j < i; j++) PetscCall(VecAXPY(gl->Z, -shift * h * a[i * s + j], Ydot[j])); 917 /* Note: Z is used within function evaluation, Ydot = Z + shift*Y */ 918 919 /* Compute an estimate of Y to start Newton iteration */ 920 if (gl->extrapolate) { 921 if (i == 0) { 922 /* Linear extrapolation on the first stage */ 923 PetscCall(VecWAXPY(Y, c[i] * h, X[1], X[0])); 924 } else { 925 /* Linear extrapolation from the last stage */ 926 PetscCall(VecAXPY(Y, (c[i] - c[i - 1]) * h, Ydot[i - 1])); 927 } 928 } else if (i == 0) { /* Directly use solution from the last step, otherwise reuse the last stage (do nothing) */ 929 PetscCall(VecCopy(X[0], Y)); 930 } 931 932 /* Solve this stage (Ydot[i] is computed during function evaluation) */ 933 PetscCall(SNESSolve(ts->snes, NULL, Y)); 934 PetscCall(SNESGetIterationNumber(ts->snes, &its)); 935 PetscCall(SNESGetLinearSolveIterations(ts->snes, &lits)); 936 PetscCall(SNESGetConvergedReason(ts->snes, &snesreason)); 937 ts->snes_its += its; 938 ts->ksp_its += lits; 939 if (snesreason < 0 && ts->max_snes_failures != PETSC_UNLIMITED && ++ts->num_snes_failures >= ts->max_snes_failures) { 940 ts->reason = TS_DIVERGED_NONLINEAR_SOLVE; 941 PetscCall(PetscInfo(ts, "Step=%" PetscInt_FMT ", nonlinear solve failures %" PetscInt_FMT " greater than current TS allowed, stopping solve\n", ts->steps, ts->num_snes_failures)); 942 PetscFunctionReturn(PETSC_SUCCESS); 943 } 944 } 945 946 gl->stage_time = ts->ptime + ts->time_step; 947 948 PetscCall((*gl->EstimateHigherMoments)(scheme, h, Ydot, gl->X, gl->himom)); 949 /* hmnorm[i] = h^{p+i}x^{(p+i)} with i=0,1,2; hmnorm[3] = h^{p+2}(dx'/dx) x^{(p+1)} */ 950 for (i = 0; i < 3; i++) PetscCall(TSGLLEVecNormWRMS(ts, gl->himom[i], &hmnorm[i + 1])); 951 enorm[0] = PetscRealPart(scheme->alpha[0]) * hmnorm[1]; 952 enorm[1] = PetscRealPart(scheme->beta[0]) * hmnorm[2]; 953 enorm[2] = PetscRealPart(scheme->gamma[0]) * hmnorm[3]; 954 PetscCall((*gl->Accept)(ts, ts->max_time - gl->stage_time, h, enorm, &accept)); 955 if (accept) goto accepted; 956 rejections++; 957 PetscCall(PetscInfo(ts, "Step %" PetscInt_FMT " (t=%g) not accepted, rejections=%" PetscInt_FMT "\n", k, (double)gl->stage_time, rejections)); 958 if (rejections > gl->max_step_rejections) break; 959 /* 960 There are lots of reasons why a step might be rejected, including solvers not converging and other factors that 961 TSGLLEChooseNextScheme does not support. Additionally, the error estimates may be very screwed up, so I'm not 962 convinced that it's safe to just compute a new error estimate using the same interface as the current adaptor 963 (the adaptor interface probably has to change). Here we make an arbitrary and naive choice. This assumes that 964 steps were written in Nordsieck form. The "correct" method would be to re-complete the previous time step with 965 the correct "next" step size. It is unclear to me whether the present ad-hoc method of rescaling X is stable. 966 */ 967 h *= 0.5; 968 for (i = 1; i < scheme->r; i++) PetscCall(VecScale(X[i], PetscPowRealInt(0.5, i))); 969 } 970 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_CONV_FAILED, "Time step %" PetscInt_FMT " (t=%g) not accepted after %" PetscInt_FMT " failures", k, (double)gl->stage_time, rejections); 971 972 accepted: 973 /* This term is not error, but it *would* be the leading term for a lower order method */ 974 PetscCall(TSGLLEVecNormWRMS(ts, gl->X[scheme->r - 1], &hmnorm[0])); 975 /* Correct scaling so that these are equivalent to norms of the Nordsieck vectors */ 976 977 PetscCall(PetscInfo(ts, "Last moment norm %10.2e, estimated error norms %10.2e %10.2e %10.2e\n", (double)hmnorm[0], (double)enorm[0], (double)enorm[1], (double)enorm[2])); 978 if (!final_step) { 979 PetscCall(TSGLLEChooseNextScheme(ts, h, hmnorm, &next_scheme, &next_h, &final_step)); 980 } else { 981 /* Dummy values to complete the current step in a consistent manner */ 982 next_scheme = gl->current_scheme; 983 next_h = h; 984 finish = PETSC_TRUE; 985 } 986 987 X = gl->Xold; 988 gl->Xold = gl->X; 989 gl->X = X; 990 PetscCall((*gl->CompleteStep)(scheme, h, gl->schemes[next_scheme], next_h, Ydot, gl->Xold, gl->X)); 991 992 PetscCall(TSGLLEUpdateWRMS(ts)); 993 994 /* Post the solution for the user, we could avoid this copy with a small bit of cleverness */ 995 PetscCall(VecCopy(gl->X[0], ts->vec_sol)); 996 ts->ptime += h; 997 ts->steps++; 998 999 PetscCall(TSPostEvaluate(ts)); 1000 PetscCall(TSPostStep(ts)); 1001 PetscCall(TSMonitor(ts, ts->steps, ts->ptime, ts->vec_sol)); 1002 1003 gl->current_scheme = next_scheme; 1004 ts->time_step = next_h; 1005 } 1006 PetscFunctionReturn(PETSC_SUCCESS); 1007 } 1008 1009 /*------------------------------------------------------------*/ 1010 1011 static PetscErrorCode TSReset_GLLE(TS ts) 1012 { 1013 TS_GLLE *gl = (TS_GLLE *)ts->data; 1014 PetscInt max_r, max_s; 1015 1016 PetscFunctionBegin; 1017 if (gl->setupcalled) { 1018 PetscCall(TSGLLEGetMaxSizes(ts, &max_r, &max_s)); 1019 PetscCall(VecDestroyVecs(max_r, &gl->Xold)); 1020 PetscCall(VecDestroyVecs(max_r, &gl->X)); 1021 PetscCall(VecDestroyVecs(max_s, &gl->Ydot)); 1022 PetscCall(VecDestroyVecs(3, &gl->himom)); 1023 PetscCall(VecDestroy(&gl->W)); 1024 PetscCall(VecDestroy(&gl->Y)); 1025 PetscCall(VecDestroy(&gl->Z)); 1026 } 1027 gl->setupcalled = PETSC_FALSE; 1028 PetscFunctionReturn(PETSC_SUCCESS); 1029 } 1030 1031 static PetscErrorCode TSDestroy_GLLE(TS ts) 1032 { 1033 TS_GLLE *gl = (TS_GLLE *)ts->data; 1034 1035 PetscFunctionBegin; 1036 PetscCall(TSReset_GLLE(ts)); 1037 if (ts->dm) { 1038 PetscCall(DMCoarsenHookRemove(ts->dm, DMCoarsenHook_TSGLLE, DMRestrictHook_TSGLLE, ts)); 1039 PetscCall(DMSubDomainHookRemove(ts->dm, DMSubDomainHook_TSGLLE, DMSubDomainRestrictHook_TSGLLE, ts)); 1040 } 1041 if (gl->adapt) PetscCall(TSGLLEAdaptDestroy(&gl->adapt)); 1042 if (gl->Destroy) PetscCall((*gl->Destroy)(gl)); 1043 PetscCall(PetscFree(ts->data)); 1044 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSGLLESetType_C", NULL)); 1045 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSGLLESetAcceptType_C", NULL)); 1046 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSGLLEGetAdapt_C", NULL)); 1047 PetscFunctionReturn(PETSC_SUCCESS); 1048 } 1049 1050 /* 1051 This defines the nonlinear equation that is to be solved with SNES 1052 g(x) = f(t,x,z+shift*x) = 0 1053 */ 1054 static PetscErrorCode SNESTSFormFunction_GLLE(SNES snes, Vec x, Vec f, TS ts) 1055 { 1056 TS_GLLE *gl = (TS_GLLE *)ts->data; 1057 Vec Z, Ydot; 1058 DM dm, dmsave; 1059 1060 PetscFunctionBegin; 1061 PetscCall(SNESGetDM(snes, &dm)); 1062 PetscCall(TSGLLEGetVecs(ts, dm, &Z, &Ydot)); 1063 PetscCall(VecWAXPY(Ydot, gl->scoeff / ts->time_step, x, Z)); 1064 dmsave = ts->dm; 1065 ts->dm = dm; 1066 PetscCall(TSComputeIFunction(ts, gl->stage_time, x, Ydot, f, PETSC_FALSE)); 1067 ts->dm = dmsave; 1068 PetscCall(TSGLLERestoreVecs(ts, dm, &Z, &Ydot)); 1069 PetscFunctionReturn(PETSC_SUCCESS); 1070 } 1071 1072 static PetscErrorCode SNESTSFormJacobian_GLLE(SNES snes, Vec x, Mat A, Mat B, TS ts) 1073 { 1074 TS_GLLE *gl = (TS_GLLE *)ts->data; 1075 Vec Z, Ydot; 1076 DM dm, dmsave; 1077 1078 PetscFunctionBegin; 1079 PetscCall(SNESGetDM(snes, &dm)); 1080 PetscCall(TSGLLEGetVecs(ts, dm, &Z, &Ydot)); 1081 dmsave = ts->dm; 1082 ts->dm = dm; 1083 /* gl->Xdot will have already been computed in SNESTSFormFunction_GLLE */ 1084 PetscCall(TSComputeIJacobian(ts, gl->stage_time, x, gl->Ydot[gl->stage], gl->scoeff / ts->time_step, A, B, PETSC_FALSE)); 1085 ts->dm = dmsave; 1086 PetscCall(TSGLLERestoreVecs(ts, dm, &Z, &Ydot)); 1087 PetscFunctionReturn(PETSC_SUCCESS); 1088 } 1089 1090 static PetscErrorCode TSSetUp_GLLE(TS ts) 1091 { 1092 TS_GLLE *gl = (TS_GLLE *)ts->data; 1093 PetscInt max_r, max_s; 1094 DM dm; 1095 1096 PetscFunctionBegin; 1097 gl->setupcalled = PETSC_TRUE; 1098 PetscCall(TSGLLEGetMaxSizes(ts, &max_r, &max_s)); 1099 PetscCall(VecDuplicateVecs(ts->vec_sol, max_r, &gl->X)); 1100 PetscCall(VecDuplicateVecs(ts->vec_sol, max_r, &gl->Xold)); 1101 PetscCall(VecDuplicateVecs(ts->vec_sol, max_s, &gl->Ydot)); 1102 PetscCall(VecDuplicateVecs(ts->vec_sol, 3, &gl->himom)); 1103 PetscCall(VecDuplicate(ts->vec_sol, &gl->W)); 1104 PetscCall(VecDuplicate(ts->vec_sol, &gl->Y)); 1105 PetscCall(VecDuplicate(ts->vec_sol, &gl->Z)); 1106 1107 /* Default acceptance tests and adaptivity */ 1108 if (!gl->Accept) PetscCall(TSGLLESetAcceptType(ts, TSGLLEACCEPT_ALWAYS)); 1109 if (!gl->adapt) PetscCall(TSGLLEGetAdapt(ts, &gl->adapt)); 1110 1111 if (gl->current_scheme < 0) { 1112 PetscInt i; 1113 for (i = 0;; i++) { 1114 if (gl->schemes[i]->p == gl->start_order) break; 1115 PetscCheck(i + 1 != gl->nschemes, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No schemes available with requested start order %" PetscInt_FMT, i); 1116 } 1117 gl->current_scheme = i; 1118 } 1119 PetscCall(TSGetDM(ts, &dm)); 1120 PetscCall(DMCoarsenHookAdd(dm, DMCoarsenHook_TSGLLE, DMRestrictHook_TSGLLE, ts)); 1121 PetscCall(DMSubDomainHookAdd(dm, DMSubDomainHook_TSGLLE, DMSubDomainRestrictHook_TSGLLE, ts)); 1122 PetscFunctionReturn(PETSC_SUCCESS); 1123 } 1124 /*------------------------------------------------------------*/ 1125 1126 static PetscErrorCode TSSetFromOptions_GLLE(TS ts, PetscOptionItems PetscOptionsObject) 1127 { 1128 TS_GLLE *gl = (TS_GLLE *)ts->data; 1129 char tname[256] = TSGLLE_IRKS, completef[256] = "rescale-and-modify"; 1130 1131 PetscFunctionBegin; 1132 PetscOptionsHeadBegin(PetscOptionsObject, "General Linear ODE solver options"); 1133 { 1134 PetscBool flg; 1135 PetscCall(PetscOptionsFList("-ts_gl_type", "Type of GL method", "TSGLLESetType", TSGLLEList, gl->type_name[0] ? gl->type_name : tname, tname, sizeof(tname), &flg)); 1136 if (flg || !gl->type_name[0]) PetscCall(TSGLLESetType(ts, tname)); 1137 PetscCall(PetscOptionsInt("-ts_gl_max_step_rejections", "Maximum number of times to attempt a step", "None", gl->max_step_rejections, &gl->max_step_rejections, NULL)); 1138 PetscCall(PetscOptionsInt("-ts_gl_max_order", "Maximum order to try", "TSGLLESetMaxOrder", gl->max_order, &gl->max_order, NULL)); 1139 PetscCall(PetscOptionsInt("-ts_gl_min_order", "Minimum order to try", "TSGLLESetMinOrder", gl->min_order, &gl->min_order, NULL)); 1140 PetscCall(PetscOptionsInt("-ts_gl_start_order", "Initial order to try", "TSGLLESetMinOrder", gl->start_order, &gl->start_order, NULL)); 1141 PetscCall(PetscOptionsEnum("-ts_gl_error_direction", "Which direction to look when estimating error", "TSGLLESetErrorDirection", TSGLLEErrorDirections, (PetscEnum)gl->error_direction, (PetscEnum *)&gl->error_direction, NULL)); 1142 PetscCall(PetscOptionsBool("-ts_gl_extrapolate", "Extrapolate stage solution from previous solution (sometimes unstable)", "TSGLLESetExtrapolate", gl->extrapolate, &gl->extrapolate, NULL)); 1143 PetscCall(PetscOptionsReal("-ts_gl_atol", "Absolute tolerance", "TSGLLESetTolerances", gl->wrms_atol, &gl->wrms_atol, NULL)); 1144 PetscCall(PetscOptionsReal("-ts_gl_rtol", "Relative tolerance", "TSGLLESetTolerances", gl->wrms_rtol, &gl->wrms_rtol, NULL)); 1145 PetscCall(PetscOptionsString("-ts_gl_complete", "Method to use for completing the step", "none", completef, completef, sizeof(completef), &flg)); 1146 if (flg) { 1147 PetscBool match1, match2; 1148 PetscCall(PetscStrcmp(completef, "rescale", &match1)); 1149 PetscCall(PetscStrcmp(completef, "rescale-and-modify", &match2)); 1150 if (match1) gl->CompleteStep = TSGLLECompleteStep_Rescale; 1151 else if (match2) gl->CompleteStep = TSGLLECompleteStep_RescaleAndModify; 1152 else SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_UNKNOWN_TYPE, "%s", completef); 1153 } 1154 { 1155 char type[256] = TSGLLEACCEPT_ALWAYS; 1156 PetscCall(PetscOptionsFList("-ts_gl_accept_type", "Method to use for determining whether to accept a step", "TSGLLESetAcceptType", TSGLLEAcceptList, gl->accept_name[0] ? gl->accept_name : type, type, sizeof(type), &flg)); 1157 if (flg || !gl->accept_name[0]) PetscCall(TSGLLESetAcceptType(ts, type)); 1158 } 1159 { 1160 TSGLLEAdapt adapt; 1161 PetscCall(TSGLLEGetAdapt(ts, &adapt)); 1162 PetscCall(TSGLLEAdaptSetFromOptions(adapt, PetscOptionsObject)); 1163 } 1164 } 1165 PetscOptionsHeadEnd(); 1166 PetscFunctionReturn(PETSC_SUCCESS); 1167 } 1168 1169 static PetscErrorCode TSView_GLLE(TS ts, PetscViewer viewer) 1170 { 1171 TS_GLLE *gl = (TS_GLLE *)ts->data; 1172 PetscInt i; 1173 PetscBool iascii, details; 1174 1175 PetscFunctionBegin; 1176 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii)); 1177 if (iascii) { 1178 PetscCall(PetscViewerASCIIPrintf(viewer, " min order %" PetscInt_FMT ", max order %" PetscInt_FMT ", current order %" PetscInt_FMT "\n", gl->min_order, gl->max_order, gl->schemes[gl->current_scheme]->p)); 1179 PetscCall(PetscViewerASCIIPrintf(viewer, " Error estimation: %s\n", TSGLLEErrorDirections[gl->error_direction])); 1180 PetscCall(PetscViewerASCIIPrintf(viewer, " Extrapolation: %s\n", gl->extrapolate ? "yes" : "no")); 1181 PetscCall(PetscViewerASCIIPrintf(viewer, " Acceptance test: %s\n", gl->accept_name[0] ? gl->accept_name : "(not yet set)")); 1182 PetscCall(PetscViewerASCIIPushTab(viewer)); 1183 PetscCall(TSGLLEAdaptView(gl->adapt, viewer)); 1184 PetscCall(PetscViewerASCIIPopTab(viewer)); 1185 PetscCall(PetscViewerASCIIPrintf(viewer, " type: %s\n", gl->type_name[0] ? gl->type_name : "(not yet set)")); 1186 PetscCall(PetscViewerASCIIPrintf(viewer, "Schemes within family (%" PetscInt_FMT "):\n", gl->nschemes)); 1187 details = PETSC_FALSE; 1188 PetscCall(PetscOptionsGetBool(((PetscObject)ts)->options, ((PetscObject)ts)->prefix, "-ts_gl_view_detailed", &details, NULL)); 1189 PetscCall(PetscViewerASCIIPushTab(viewer)); 1190 for (i = 0; i < gl->nschemes; i++) PetscCall(TSGLLESchemeView(gl->schemes[i], details, viewer)); 1191 if (gl->View) PetscCall((*gl->View)(gl, viewer)); 1192 PetscCall(PetscViewerASCIIPopTab(viewer)); 1193 } 1194 PetscFunctionReturn(PETSC_SUCCESS); 1195 } 1196 1197 /*@C 1198 TSGLLERegister - adds a `TSGLLE` implementation 1199 1200 Not Collective, No Fortran Support 1201 1202 Input Parameters: 1203 + sname - name of user-defined general linear scheme 1204 - function - routine to create method context 1205 1206 Level: advanced 1207 1208 Note: 1209 `TSGLLERegister()` may be called multiple times to add several user-defined families. 1210 1211 Example Usage: 1212 .vb 1213 TSGLLERegister("my_scheme", MySchemeCreate); 1214 .ve 1215 1216 Then, your scheme can be chosen with the procedural interface via 1217 $ TSGLLESetType(ts, "my_scheme") 1218 or at runtime via the option 1219 $ -ts_gl_type my_scheme 1220 1221 .seealso: [](ch_ts), `TSGLLE`, `TSGLLEType`, `TSGLLERegisterAll()` 1222 @*/ 1223 PetscErrorCode TSGLLERegister(const char sname[], PetscErrorCode (*function)(TS)) 1224 { 1225 PetscFunctionBegin; 1226 PetscCall(TSGLLEInitializePackage()); 1227 PetscCall(PetscFunctionListAdd(&TSGLLEList, sname, function)); 1228 PetscFunctionReturn(PETSC_SUCCESS); 1229 } 1230 1231 /*@C 1232 TSGLLEAcceptRegister - adds a `TSGLLE` acceptance scheme 1233 1234 Not Collective 1235 1236 Input Parameters: 1237 + sname - name of user-defined acceptance scheme 1238 - function - routine to create method context, see `TSGLLEAcceptFn` for the calling sequence 1239 1240 Level: advanced 1241 1242 Note: 1243 `TSGLLEAcceptRegister()` may be called multiple times to add several user-defined families. 1244 1245 Example Usage: 1246 .vb 1247 TSGLLEAcceptRegister("my_scheme", MySchemeCreate); 1248 .ve 1249 1250 Then, your scheme can be chosen with the procedural interface via 1251 .vb 1252 TSGLLESetAcceptType(ts, "my_scheme") 1253 .ve 1254 or at runtime via the option `-ts_gl_accept_type my_scheme` 1255 1256 .seealso: [](ch_ts), `TSGLLE`, `TSGLLEType`, `TSGLLERegisterAll()`, `TSGLLEAcceptFn` 1257 @*/ 1258 PetscErrorCode TSGLLEAcceptRegister(const char sname[], TSGLLEAcceptFn *function) 1259 { 1260 PetscFunctionBegin; 1261 PetscCall(PetscFunctionListAdd(&TSGLLEAcceptList, sname, function)); 1262 PetscFunctionReturn(PETSC_SUCCESS); 1263 } 1264 1265 /*@C 1266 TSGLLERegisterAll - Registers all of the general linear methods in `TSGLLE` 1267 1268 Not Collective 1269 1270 Level: advanced 1271 1272 .seealso: [](ch_ts), `TSGLLE`, `TSGLLERegisterDestroy()` 1273 @*/ 1274 PetscErrorCode TSGLLERegisterAll(void) 1275 { 1276 PetscFunctionBegin; 1277 if (TSGLLERegisterAllCalled) PetscFunctionReturn(PETSC_SUCCESS); 1278 TSGLLERegisterAllCalled = PETSC_TRUE; 1279 1280 PetscCall(TSGLLERegister(TSGLLE_IRKS, TSGLLECreate_IRKS)); 1281 PetscCall(TSGLLEAcceptRegister(TSGLLEACCEPT_ALWAYS, TSGLLEAccept_Always)); 1282 PetscFunctionReturn(PETSC_SUCCESS); 1283 } 1284 1285 /*@C 1286 TSGLLEInitializePackage - This function initializes everything in the `TSGLLE` package. It is called 1287 from `TSInitializePackage()`. 1288 1289 Level: developer 1290 1291 .seealso: [](ch_ts), `PetscInitialize()`, `TSInitializePackage()`, `TSGLLEFinalizePackage()` 1292 @*/ 1293 PetscErrorCode TSGLLEInitializePackage(void) 1294 { 1295 PetscFunctionBegin; 1296 if (TSGLLEPackageInitialized) PetscFunctionReturn(PETSC_SUCCESS); 1297 TSGLLEPackageInitialized = PETSC_TRUE; 1298 PetscCall(TSGLLERegisterAll()); 1299 PetscCall(PetscRegisterFinalize(TSGLLEFinalizePackage)); 1300 PetscFunctionReturn(PETSC_SUCCESS); 1301 } 1302 1303 /*@C 1304 TSGLLEFinalizePackage - This function destroys everything in the `TSGLLE` package. It is 1305 called from `PetscFinalize()`. 1306 1307 Level: developer 1308 1309 .seealso: [](ch_ts), `PetscFinalize()`, `TSGLLEInitializePackage()`, `TSInitializePackage()` 1310 @*/ 1311 PetscErrorCode TSGLLEFinalizePackage(void) 1312 { 1313 PetscFunctionBegin; 1314 PetscCall(PetscFunctionListDestroy(&TSGLLEList)); 1315 PetscCall(PetscFunctionListDestroy(&TSGLLEAcceptList)); 1316 TSGLLEPackageInitialized = PETSC_FALSE; 1317 TSGLLERegisterAllCalled = PETSC_FALSE; 1318 PetscFunctionReturn(PETSC_SUCCESS); 1319 } 1320 1321 /* ------------------------------------------------------------ */ 1322 /*MC 1323 TSGLLE - DAE solver using implicit General Linear methods {cite}`butcher_2007` {cite}`butcher2016numerical` 1324 1325 Options Database Keys: 1326 + -ts_gl_type <type> - the class of general linear method (irks) 1327 . -ts_gl_rtol <tol> - relative error 1328 . -ts_gl_atol <tol> - absolute error 1329 . -ts_gl_min_order <p> - minimum order method to consider (default=1) 1330 . -ts_gl_max_order <p> - maximum order method to consider (default=3) 1331 . -ts_gl_start_order <p> - order of starting method (default=1) 1332 . -ts_gl_complete <method> - method to use for completing the step (rescale-and-modify or rescale) 1333 - -ts_adapt_type <method> - adaptive controller to use (none step both) 1334 1335 Level: beginner 1336 1337 Notes: 1338 These methods contain Runge-Kutta and multistep schemes as special cases. These special cases 1339 have some fundamental limitations. For example, diagonally implicit Runge-Kutta cannot have 1340 stage order greater than 1 which limits their applicability to very stiff systems. 1341 Meanwhile, multistep methods cannot be A-stable for order greater than 2 and BDF are not 1342 0-stable for order greater than 6. GL methods can be A- and L-stable with arbitrarily high 1343 stage order and reliable error estimates for both 1 and 2 orders higher to facilitate 1344 adaptive step sizes and adaptive order schemes. All this is possible while preserving a 1345 singly diagonally implicit structure. 1346 1347 This integrator can be applied to DAE. 1348 1349 Diagonally implicit general linear (DIGL) methods are a generalization of diagonally implicit 1350 Runge-Kutta (DIRK). They are represented by the tableau 1351 1352 .vb 1353 A | U 1354 ------- 1355 B | V 1356 .ve 1357 1358 combined with a vector c of abscissa. "Diagonally implicit" means that $A$ is lower 1359 triangular. A step of the general method reads 1360 1361 $$ 1362 \begin{align*} 1363 [ Y ] = [A U] [ Y' ] \\ 1364 [X^k] = [B V] [X^{k-1}] 1365 \end{align*} 1366 $$ 1367 1368 where Y is the multivector of stage values, $Y'$ is the multivector of stage derivatives, $X^k$ 1369 is the Nordsieck vector of the solution at step $k$. The Nordsieck vector consists of the first 1370 $r$ moments of the solution, given by 1371 1372 $$ 1373 X = [x_0,x_1,...,x_{r-1}] = [x, h x', h^2 x'', ..., h^{r-1} x^{(r-1)} ] 1374 $$ 1375 1376 If $A$ is lower triangular, we can solve the stages $(Y, Y')$ sequentially 1377 1378 $$ 1379 y_i = h \sum_{j=0}^{s-1} (a_{ij} y'_j) + \sum_{j=0}^{r-1} u_{ij} x_j, \, \, i=0,...,{s-1} 1380 $$ 1381 1382 and then construct the pieces to carry to the next step 1383 1384 $$ 1385 xx_i = h \sum_{j=0}^{s-1} b_{ij} y'_j + \sum_{j=0}^{r-1} v_{ij} x_j, \, \, i=0,...,{r-1} 1386 $$ 1387 1388 Note that when the equations are cast in implicit form, we are using the stage equation to 1389 define $y'_i$ in terms of $y_i$ and known stuff ($y_j$ for $j<i$ and $x_j$ for all $j$). 1390 1391 Error estimation 1392 1393 At present, the most attractive GL methods for stiff problems are singly diagonally implicit 1394 schemes which posses Inherent Runge-Kutta Stability (`TSIRKS`). These methods have $r=s$, the 1395 number of items passed between steps is equal to the number of stages. The order and 1396 stage-order are one less than the number of stages. We use the error estimates in the 2007 1397 paper which provide the following estimates 1398 1399 $$ 1400 \begin{align*} 1401 h^{p+1} X^{(p+1)} = \phi_0^T Y' + [0 \psi_0^T] Xold \\ 1402 h^{p+2} X^{(p+2)} = \phi_1^T Y' + [0 \psi_1^T] Xold \\ 1403 h^{p+2} (dx'/dx) X^{(p+1)} = \phi_2^T Y' + [0 \psi_2^T] Xold 1404 \end{align*} 1405 $$ 1406 1407 These estimates are accurate to $ O(h^{p+3})$. 1408 1409 Changing the step size 1410 1411 Uses the generalized "rescale and modify" scheme, see equation (4.5) of {cite}`butcher_2007`. 1412 1413 .seealso: [](ch_ts), `TSCreate()`, `TS`, `TSSetType()`, `TSType` 1414 M*/ 1415 PETSC_EXTERN PetscErrorCode TSCreate_GLLE(TS ts) 1416 { 1417 TS_GLLE *gl; 1418 1419 PetscFunctionBegin; 1420 PetscCall(TSGLLEInitializePackage()); 1421 1422 PetscCall(PetscNew(&gl)); 1423 ts->data = (void *)gl; 1424 1425 ts->ops->reset = TSReset_GLLE; 1426 ts->ops->destroy = TSDestroy_GLLE; 1427 ts->ops->view = TSView_GLLE; 1428 ts->ops->setup = TSSetUp_GLLE; 1429 ts->ops->solve = TSSolve_GLLE; 1430 ts->ops->setfromoptions = TSSetFromOptions_GLLE; 1431 ts->ops->snesfunction = SNESTSFormFunction_GLLE; 1432 ts->ops->snesjacobian = SNESTSFormJacobian_GLLE; 1433 1434 ts->usessnes = PETSC_TRUE; 1435 1436 gl->max_step_rejections = 1; 1437 gl->min_order = 1; 1438 gl->max_order = 3; 1439 gl->start_order = 1; 1440 gl->current_scheme = -1; 1441 gl->extrapolate = PETSC_FALSE; 1442 1443 gl->wrms_atol = 1e-8; 1444 gl->wrms_rtol = 1e-5; 1445 1446 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSGLLESetType_C", &TSGLLESetType_GLLE)); 1447 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSGLLESetAcceptType_C", &TSGLLESetAcceptType_GLLE)); 1448 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSGLLEGetAdapt_C", &TSGLLEGetAdapt_GLLE)); 1449 PetscFunctionReturn(PETSC_SUCCESS); 1450 } 1451