1 2 #include <../src/ts/impls/implicit/glle/glle.h> /*I "petscts.h" I*/ 3 #include <petscdm.h> 4 #include <petscblaslapack.h> 5 6 static const char *TSGLLEErrorDirections[] = {"FORWARD", "BACKWARD", "TSGLLEErrorDirection", "TSGLLEERROR_", NULL}; 7 static PetscFunctionList TSGLLEList; 8 static PetscFunctionList TSGLLEAcceptList; 9 static PetscBool TSGLLEPackageInitialized; 10 static PetscBool TSGLLERegisterAllCalled; 11 12 /* This function is pure */ 13 static PetscScalar Factorial(PetscInt n) 14 { 15 PetscInt i; 16 if (n < 12) { /* Can compute with 32-bit integers */ 17 PetscInt f = 1; 18 for (i = 2; i <= n; i++) f *= i; 19 return (PetscScalar)f; 20 } else { 21 PetscScalar f = 1.; 22 for (i = 2; i <= n; i++) f *= (PetscScalar)i; 23 return f; 24 } 25 } 26 27 /* This function is pure */ 28 static PetscScalar CPowF(PetscScalar c, PetscInt p) 29 { 30 return PetscPowRealInt(PetscRealPart(c), p) / Factorial(p); 31 } 32 33 static PetscErrorCode TSGLLEGetVecs(TS ts, DM dm, Vec *Z, Vec *Ydotstage) 34 { 35 TS_GLLE *gl = (TS_GLLE *)ts->data; 36 37 PetscFunctionBegin; 38 if (Z) { 39 if (dm && dm != ts->dm) { 40 PetscCall(DMGetNamedGlobalVector(dm, "TSGLLE_Z", Z)); 41 } else *Z = gl->Z; 42 } 43 if (Ydotstage) { 44 if (dm && dm != ts->dm) { 45 PetscCall(DMGetNamedGlobalVector(dm, "TSGLLE_Ydot", Ydotstage)); 46 } else *Ydotstage = gl->Ydot[gl->stage]; 47 } 48 PetscFunctionReturn(0); 49 } 50 51 static PetscErrorCode TSGLLERestoreVecs(TS ts, DM dm, Vec *Z, Vec *Ydotstage) 52 { 53 PetscFunctionBegin; 54 if (Z) { 55 if (dm && dm != ts->dm) PetscCall(DMRestoreNamedGlobalVector(dm, "TSGLLE_Z", Z)); 56 } 57 if (Ydotstage) { 58 if (dm && dm != ts->dm) PetscCall(DMRestoreNamedGlobalVector(dm, "TSGLLE_Ydot", Ydotstage)); 59 } 60 PetscFunctionReturn(0); 61 } 62 63 static PetscErrorCode DMCoarsenHook_TSGLLE(DM fine, DM coarse, void *ctx) 64 { 65 PetscFunctionBegin; 66 PetscFunctionReturn(0); 67 } 68 69 static PetscErrorCode DMRestrictHook_TSGLLE(DM fine, Mat restrct, Vec rscale, Mat inject, DM coarse, void *ctx) 70 { 71 TS ts = (TS)ctx; 72 Vec Ydot, Ydot_c; 73 74 PetscFunctionBegin; 75 PetscCall(TSGLLEGetVecs(ts, fine, NULL, &Ydot)); 76 PetscCall(TSGLLEGetVecs(ts, coarse, NULL, &Ydot_c)); 77 PetscCall(MatRestrict(restrct, Ydot, Ydot_c)); 78 PetscCall(VecPointwiseMult(Ydot_c, rscale, Ydot_c)); 79 PetscCall(TSGLLERestoreVecs(ts, fine, NULL, &Ydot)); 80 PetscCall(TSGLLERestoreVecs(ts, coarse, NULL, &Ydot_c)); 81 PetscFunctionReturn(0); 82 } 83 84 static PetscErrorCode DMSubDomainHook_TSGLLE(DM dm, DM subdm, void *ctx) 85 { 86 PetscFunctionBegin; 87 PetscFunctionReturn(0); 88 } 89 90 static PetscErrorCode DMSubDomainRestrictHook_TSGLLE(DM dm, VecScatter gscat, VecScatter lscat, DM subdm, void *ctx) 91 { 92 TS ts = (TS)ctx; 93 Vec Ydot, Ydot_s; 94 95 PetscFunctionBegin; 96 PetscCall(TSGLLEGetVecs(ts, dm, NULL, &Ydot)); 97 PetscCall(TSGLLEGetVecs(ts, subdm, NULL, &Ydot_s)); 98 99 PetscCall(VecScatterBegin(gscat, Ydot, Ydot_s, INSERT_VALUES, SCATTER_FORWARD)); 100 PetscCall(VecScatterEnd(gscat, Ydot, Ydot_s, INSERT_VALUES, SCATTER_FORWARD)); 101 102 PetscCall(TSGLLERestoreVecs(ts, dm, NULL, &Ydot)); 103 PetscCall(TSGLLERestoreVecs(ts, subdm, NULL, &Ydot_s)); 104 PetscFunctionReturn(0); 105 } 106 107 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) 108 { 109 TSGLLEScheme scheme; 110 PetscInt j; 111 112 PetscFunctionBegin; 113 PetscCheck(p >= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Scheme order must be positive"); 114 PetscCheck(r >= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "At least one item must be carried between steps"); 115 PetscCheck(s >= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "At least one stage is required"); 116 PetscValidPointer(inscheme, 10); 117 *inscheme = NULL; 118 PetscCall(PetscNew(&scheme)); 119 scheme->p = p; 120 scheme->q = q; 121 scheme->r = r; 122 scheme->s = s; 123 124 PetscCall(PetscMalloc5(s, &scheme->c, s * s, &scheme->a, r * s, &scheme->b, r * s, &scheme->u, r * r, &scheme->v)); 125 PetscCall(PetscArraycpy(scheme->c, c, s)); 126 for (j = 0; j < s * s; j++) scheme->a[j] = (PetscAbsScalar(a[j]) < 1e-12) ? 0 : a[j]; 127 for (j = 0; j < r * s; j++) scheme->b[j] = (PetscAbsScalar(b[j]) < 1e-12) ? 0 : b[j]; 128 for (j = 0; j < s * r; j++) scheme->u[j] = (PetscAbsScalar(u[j]) < 1e-12) ? 0 : u[j]; 129 for (j = 0; j < r * r; j++) scheme->v[j] = (PetscAbsScalar(v[j]) < 1e-12) ? 0 : v[j]; 130 131 PetscCall(PetscMalloc6(r, &scheme->alpha, r, &scheme->beta, r, &scheme->gamma, 3 * s, &scheme->phi, 3 * r, &scheme->psi, r, &scheme->stage_error)); 132 { 133 PetscInt i, j, k, ss = s + 2; 134 PetscBLASInt m, n, one = 1, *ipiv, lwork = 4 * ((s + 3) * 3 + 3), info, ldb; 135 PetscReal rcond, *sing, *workreal; 136 PetscScalar *ImV, *H, *bmat, *workscalar, *c = scheme->c, *a = scheme->a, *b = scheme->b, *u = scheme->u, *v = scheme->v; 137 PetscBLASInt rank; 138 PetscCall(PetscMalloc7(PetscSqr(r), &ImV, 3 * s, &H, 3 * ss, &bmat, lwork, &workscalar, 5 * (3 + r), &workreal, r + s, &sing, r + s, &ipiv)); 139 140 /* column-major input */ 141 for (i = 0; i < r - 1; i++) { 142 for (j = 0; j < r - 1; j++) ImV[i + j * r] = 1.0 * (i == j) - v[(i + 1) * r + j + 1]; 143 } 144 /* Build right hand side for alpha (tp - glm.B(2:end,:)*(glm.c.^(p)./factorial(p))) */ 145 for (i = 1; i < r; i++) { 146 scheme->alpha[i] = 1. / Factorial(p + 1 - i); 147 for (j = 0; j < s; j++) scheme->alpha[i] -= b[i * s + j] * CPowF(c[j], p); 148 } 149 PetscCall(PetscBLASIntCast(r - 1, &m)); 150 PetscCall(PetscBLASIntCast(r, &n)); 151 PetscCallBLAS("LAPACKgesv", LAPACKgesv_(&m, &one, ImV, &n, ipiv, scheme->alpha + 1, &n, &info)); 152 PetscCheck(info >= 0, PETSC_COMM_SELF, PETSC_ERR_LIB, "Bad argument to GESV"); 153 PetscCheck(info <= 0, PETSC_COMM_SELF, PETSC_ERR_MAT_LU_ZRPVT, "Bad LU factorization"); 154 155 /* Build right hand side for beta (tp1 - glm.B(2:end,:)*(glm.c.^(p+1)./factorial(p+1)) - e.alpha) */ 156 for (i = 1; i < r; i++) { 157 scheme->beta[i] = 1. / Factorial(p + 2 - i) - scheme->alpha[i]; 158 for (j = 0; j < s; j++) scheme->beta[i] -= b[i * s + j] * CPowF(c[j], p + 1); 159 } 160 PetscCallBLAS("LAPACKgetrs", LAPACKgetrs_("No transpose", &m, &one, ImV, &n, ipiv, scheme->beta + 1, &n, &info)); 161 PetscCheck(info >= 0, PETSC_COMM_SELF, PETSC_ERR_LIB, "Bad argument to GETRS"); 162 PetscCheck(info <= 0, PETSC_COMM_SELF, PETSC_ERR_LIB, "Should not happen"); 163 164 /* Build stage_error vector 165 xi = glm.c.^(p+1)/factorial(p+1) - glm.A*glm.c.^p/factorial(p) + glm.U(:,2:end)*e.alpha; 166 */ 167 for (i = 0; i < s; i++) { 168 scheme->stage_error[i] = CPowF(c[i], p + 1); 169 for (j = 0; j < s; j++) scheme->stage_error[i] -= a[i * s + j] * CPowF(c[j], p); 170 for (j = 1; j < r; j++) scheme->stage_error[i] += u[i * r + j] * scheme->alpha[j]; 171 } 172 173 /* alpha[0] (epsilon in B,J,W 2007) 174 epsilon = 1/factorial(p+1) - B(1,:)*c.^p/factorial(p) + V(1,2:end)*e.alpha; 175 */ 176 scheme->alpha[0] = 1. / Factorial(p + 1); 177 for (j = 0; j < s; j++) scheme->alpha[0] -= b[0 * s + j] * CPowF(c[j], p); 178 for (j = 1; j < r; j++) scheme->alpha[0] += v[0 * r + j] * scheme->alpha[j]; 179 180 /* right hand side for gamma (glm.B(2:end,:)*e.xi - e.epsilon*eye(s-1,1)) */ 181 for (i = 1; i < r; i++) { 182 scheme->gamma[i] = (i == 1 ? -1. : 0) * scheme->alpha[0]; 183 for (j = 0; j < s; j++) scheme->gamma[i] += b[i * s + j] * scheme->stage_error[j]; 184 } 185 PetscCallBLAS("LAPACKgetrs", LAPACKgetrs_("No transpose", &m, &one, ImV, &n, ipiv, scheme->gamma + 1, &n, &info)); 186 PetscCheck(info >= 0, PETSC_COMM_SELF, PETSC_ERR_LIB, "Bad argument to GETRS"); 187 PetscCheck(info <= 0, PETSC_COMM_SELF, PETSC_ERR_LIB, "Should not happen"); 188 189 /* beta[0] (rho in B,J,W 2007) 190 e.rho = 1/factorial(p+2) - glm.B(1,:)*glm.c.^(p+1)/factorial(p+1) ... 191 + glm.V(1,2:end)*e.beta;% - e.epsilon; 192 % Note: The paper (B,J,W 2007) includes the last term in their definition 193 * */ 194 scheme->beta[0] = 1. / Factorial(p + 2); 195 for (j = 0; j < s; j++) scheme->beta[0] -= b[0 * s + j] * CPowF(c[j], p + 1); 196 for (j = 1; j < r; j++) scheme->beta[0] += v[0 * r + j] * scheme->beta[j]; 197 198 /* gamma[0] (sigma in B,J,W 2007) 199 * e.sigma = glm.B(1,:)*e.xi + glm.V(1,2:end)*e.gamma; 200 * */ 201 scheme->gamma[0] = 0.0; 202 for (j = 0; j < s; j++) scheme->gamma[0] += b[0 * s + j] * scheme->stage_error[j]; 203 for (j = 1; j < r; j++) scheme->gamma[0] += v[0 * s + j] * scheme->gamma[j]; 204 205 /* Assemble H 206 * % " PetscInt_FMT "etermine the error estimators phi 207 H = [[cpow(glm.c,p) + C*e.alpha] [cpow(glm.c,p+1) + C*e.beta] ... 208 [e.xi - C*(e.gamma + 0*e.epsilon*eye(s-1,1))]]'; 209 % Paper has formula above without the 0, but that term must be left 210 % out to satisfy the conditions they propose and to make the 211 % example schemes work 212 e.H = H; 213 e.phi = (H \ [1 0 0;1 1 0;0 0 -1])'; 214 e.psi = -e.phi*C; 215 * */ 216 for (j = 0; j < s; j++) { 217 H[0 + j * 3] = CPowF(c[j], p); 218 H[1 + j * 3] = CPowF(c[j], p + 1); 219 H[2 + j * 3] = scheme->stage_error[j]; 220 for (k = 1; k < r; k++) { 221 H[0 + j * 3] += CPowF(c[j], k - 1) * scheme->alpha[k]; 222 H[1 + j * 3] += CPowF(c[j], k - 1) * scheme->beta[k]; 223 H[2 + j * 3] -= CPowF(c[j], k - 1) * scheme->gamma[k]; 224 } 225 } 226 bmat[0 + 0 * ss] = 1.; 227 bmat[0 + 1 * ss] = 0.; 228 bmat[0 + 2 * ss] = 0.; 229 bmat[1 + 0 * ss] = 1.; 230 bmat[1 + 1 * ss] = 1.; 231 bmat[1 + 2 * ss] = 0.; 232 bmat[2 + 0 * ss] = 0.; 233 bmat[2 + 1 * ss] = 0.; 234 bmat[2 + 2 * ss] = -1.; 235 m = 3; 236 PetscCall(PetscBLASIntCast(s, &n)); 237 PetscCall(PetscBLASIntCast(ss, &ldb)); 238 rcond = 1e-12; 239 #if defined(PETSC_USE_COMPLEX) 240 /* ZGELSS( M, N, NRHS, A, LDA, B, LDB, S, RCOND, RANK, WORK, LWORK, RWORK, INFO) */ 241 PetscCallBLAS("LAPACKgelss", LAPACKgelss_(&m, &n, &m, H, &m, bmat, &ldb, sing, &rcond, &rank, workscalar, &lwork, workreal, &info)); 242 #else 243 /* DGELSS( M, N, NRHS, A, LDA, B, LDB, S, RCOND, RANK, WORK, LWORK, INFO) */ 244 PetscCallBLAS("LAPACKgelss", LAPACKgelss_(&m, &n, &m, H, &m, bmat, &ldb, sing, &rcond, &rank, workscalar, &lwork, &info)); 245 #endif 246 PetscCheck(info >= 0, PETSC_COMM_SELF, PETSC_ERR_LIB, "Bad argument to GELSS"); 247 PetscCheck(info <= 0, PETSC_COMM_SELF, PETSC_ERR_LIB, "SVD failed to converge"); 248 249 for (j = 0; j < 3; j++) { 250 for (k = 0; k < s; k++) scheme->phi[k + j * s] = bmat[k + j * ss]; 251 } 252 253 /* the other part of the error estimator, psi in B,J,W 2007 */ 254 scheme->psi[0 * r + 0] = 0.; 255 scheme->psi[1 * r + 0] = 0.; 256 scheme->psi[2 * r + 0] = 0.; 257 for (j = 1; j < r; j++) { 258 scheme->psi[0 * r + j] = 0.; 259 scheme->psi[1 * r + j] = 0.; 260 scheme->psi[2 * r + j] = 0.; 261 for (k = 0; k < s; k++) { 262 scheme->psi[0 * r + j] -= CPowF(c[k], j - 1) * scheme->phi[0 * s + k]; 263 scheme->psi[1 * r + j] -= CPowF(c[k], j - 1) * scheme->phi[1 * s + k]; 264 scheme->psi[2 * r + j] -= CPowF(c[k], j - 1) * scheme->phi[2 * s + k]; 265 } 266 } 267 PetscCall(PetscFree7(ImV, H, bmat, workscalar, workreal, sing, ipiv)); 268 } 269 /* Check which properties are satisfied */ 270 scheme->stiffly_accurate = PETSC_TRUE; 271 if (scheme->c[s - 1] != 1.) scheme->stiffly_accurate = PETSC_FALSE; 272 for (j = 0; j < s; j++) 273 if (a[(s - 1) * s + j] != b[j]) scheme->stiffly_accurate = PETSC_FALSE; 274 for (j = 0; j < r; j++) 275 if (u[(s - 1) * r + j] != v[j]) scheme->stiffly_accurate = PETSC_FALSE; 276 scheme->fsal = scheme->stiffly_accurate; /* FSAL is stronger */ 277 for (j = 0; j < s - 1; j++) 278 if (r > 1 && b[1 * s + j] != 0.) scheme->fsal = PETSC_FALSE; 279 if (b[1 * s + r - 1] != 1.) scheme->fsal = PETSC_FALSE; 280 for (j = 0; j < r; j++) 281 if (r > 1 && v[1 * r + j] != 0.) scheme->fsal = PETSC_FALSE; 282 283 *inscheme = scheme; 284 PetscFunctionReturn(0); 285 } 286 287 static PetscErrorCode TSGLLESchemeDestroy(TSGLLEScheme sc) 288 { 289 PetscFunctionBegin; 290 PetscCall(PetscFree5(sc->c, sc->a, sc->b, sc->u, sc->v)); 291 PetscCall(PetscFree6(sc->alpha, sc->beta, sc->gamma, sc->phi, sc->psi, sc->stage_error)); 292 PetscCall(PetscFree(sc)); 293 PetscFunctionReturn(0); 294 } 295 296 static PetscErrorCode TSGLLEDestroy_Default(TS_GLLE *gl) 297 { 298 PetscInt i; 299 300 PetscFunctionBegin; 301 for (i = 0; i < gl->nschemes; i++) { 302 if (gl->schemes[i]) PetscCall(TSGLLESchemeDestroy(gl->schemes[i])); 303 } 304 PetscCall(PetscFree(gl->schemes)); 305 gl->nschemes = 0; 306 PetscCall(PetscMemzero(gl->type_name, sizeof(gl->type_name))); 307 PetscFunctionReturn(0); 308 } 309 310 static PetscErrorCode TSGLLEViewTable_Private(PetscViewer viewer, PetscInt m, PetscInt n, const PetscScalar a[], const char name[]) 311 { 312 PetscBool iascii; 313 PetscInt i, j; 314 315 PetscFunctionBegin; 316 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii)); 317 if (iascii) { 318 PetscCall(PetscViewerASCIIPrintf(viewer, "%30s = [", name)); 319 for (i = 0; i < m; i++) { 320 if (i) PetscCall(PetscViewerASCIIPrintf(viewer, "%30s [", "")); 321 PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_FALSE)); 322 for (j = 0; j < n; j++) PetscCall(PetscViewerASCIIPrintf(viewer, " %12.8g", (double)PetscRealPart(a[i * n + j]))); 323 PetscCall(PetscViewerASCIIPrintf(viewer, "]\n")); 324 PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_TRUE)); 325 } 326 } 327 PetscFunctionReturn(0); 328 } 329 330 static PetscErrorCode TSGLLESchemeView(TSGLLEScheme sc, PetscBool view_details, PetscViewer viewer) 331 { 332 PetscBool iascii; 333 334 PetscFunctionBegin; 335 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii)); 336 if (iascii) { 337 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)); 338 PetscCall(PetscViewerASCIIPushTab(viewer)); 339 PetscCall(PetscViewerASCIIPrintf(viewer, "Stiffly accurate: %s, FSAL: %s\n", sc->stiffly_accurate ? "yes" : "no", sc->fsal ? "yes" : "no")); 340 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]))); 341 PetscCall(TSGLLEViewTable_Private(viewer, 1, sc->s, sc->c, "Abscissas c")); 342 if (view_details) { 343 PetscCall(TSGLLEViewTable_Private(viewer, sc->s, sc->s, sc->a, "A")); 344 PetscCall(TSGLLEViewTable_Private(viewer, sc->r, sc->s, sc->b, "B")); 345 PetscCall(TSGLLEViewTable_Private(viewer, sc->s, sc->r, sc->u, "U")); 346 PetscCall(TSGLLEViewTable_Private(viewer, sc->r, sc->r, sc->v, "V")); 347 348 PetscCall(TSGLLEViewTable_Private(viewer, 3, sc->s, sc->phi, "Error estimate phi")); 349 PetscCall(TSGLLEViewTable_Private(viewer, 3, sc->r, sc->psi, "Error estimate psi")); 350 PetscCall(TSGLLEViewTable_Private(viewer, 1, sc->r, sc->alpha, "Modify alpha")); 351 PetscCall(TSGLLEViewTable_Private(viewer, 1, sc->r, sc->beta, "Modify beta")); 352 PetscCall(TSGLLEViewTable_Private(viewer, 1, sc->r, sc->gamma, "Modify gamma")); 353 PetscCall(TSGLLEViewTable_Private(viewer, 1, sc->s, sc->stage_error, "Stage error xi")); 354 } 355 PetscCall(PetscViewerASCIIPopTab(viewer)); 356 } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Viewer type %s not supported", ((PetscObject)viewer)->type_name); 357 PetscFunctionReturn(0); 358 } 359 360 static PetscErrorCode TSGLLEEstimateHigherMoments_Default(TSGLLEScheme sc, PetscReal h, Vec Ydot[], Vec Xold[], Vec hm[]) 361 { 362 PetscInt i; 363 364 PetscFunctionBegin; 365 PetscCheck(sc->r <= 64 && sc->s <= 64, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Ridiculous number of stages or items passed between stages"); 366 /* build error vectors*/ 367 for (i = 0; i < 3; i++) { 368 PetscScalar phih[64]; 369 PetscInt j; 370 for (j = 0; j < sc->s; j++) phih[j] = sc->phi[i * sc->s + j] * h; 371 PetscCall(VecZeroEntries(hm[i])); 372 PetscCall(VecMAXPY(hm[i], sc->s, phih, Ydot)); 373 PetscCall(VecMAXPY(hm[i], sc->r, &sc->psi[i * sc->r], Xold)); 374 } 375 PetscFunctionReturn(0); 376 } 377 378 static PetscErrorCode TSGLLECompleteStep_Rescale(TSGLLEScheme sc, PetscReal h, TSGLLEScheme next_sc, PetscReal next_h, Vec Ydot[], Vec Xold[], Vec X[]) 379 { 380 PetscScalar brow[32], vrow[32]; 381 PetscInt i, j, r, s; 382 383 PetscFunctionBegin; 384 /* Build the new solution from (X,Ydot) */ 385 r = sc->r; 386 s = sc->s; 387 for (i = 0; i < r; i++) { 388 PetscCall(VecZeroEntries(X[i])); 389 for (j = 0; j < s; j++) brow[j] = h * sc->b[i * s + j]; 390 PetscCall(VecMAXPY(X[i], s, brow, Ydot)); 391 for (j = 0; j < r; j++) vrow[j] = sc->v[i * r + j]; 392 PetscCall(VecMAXPY(X[i], r, vrow, Xold)); 393 } 394 PetscFunctionReturn(0); 395 } 396 397 static PetscErrorCode TSGLLECompleteStep_RescaleAndModify(TSGLLEScheme sc, PetscReal h, TSGLLEScheme next_sc, PetscReal next_h, Vec Ydot[], Vec Xold[], Vec X[]) 398 { 399 PetscScalar brow[32], vrow[32]; 400 PetscReal ratio; 401 PetscInt i, j, p, r, s; 402 403 PetscFunctionBegin; 404 /* Build the new solution from (X,Ydot) */ 405 p = sc->p; 406 r = sc->r; 407 s = sc->s; 408 ratio = next_h / h; 409 for (i = 0; i < r; i++) { 410 PetscCall(VecZeroEntries(X[i])); 411 for (j = 0; j < s; j++) { 412 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])); 413 } 414 PetscCall(VecMAXPY(X[i], s, brow, Ydot)); 415 for (j = 0; j < r; j++) { 416 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])); 417 } 418 PetscCall(VecMAXPY(X[i], r, vrow, Xold)); 419 } 420 if (r < next_sc->r) { 421 PetscCheck(r + 1 == next_sc->r, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Cannot accommodate jump in r greater than 1"); 422 PetscCall(VecZeroEntries(X[r])); 423 for (j = 0; j < s; j++) brow[j] = h * PetscPowRealInt(ratio, p + 1) * sc->phi[0 * s + j]; 424 PetscCall(VecMAXPY(X[r], s, brow, Ydot)); 425 for (j = 0; j < r; j++) vrow[j] = PetscPowRealInt(ratio, p + 1) * sc->psi[0 * r + j]; 426 PetscCall(VecMAXPY(X[r], r, vrow, Xold)); 427 } 428 PetscFunctionReturn(0); 429 } 430 431 static PetscErrorCode TSGLLECreate_IRKS(TS ts) 432 { 433 TS_GLLE *gl = (TS_GLLE *)ts->data; 434 435 PetscFunctionBegin; 436 gl->Destroy = TSGLLEDestroy_Default; 437 gl->EstimateHigherMoments = TSGLLEEstimateHigherMoments_Default; 438 gl->CompleteStep = TSGLLECompleteStep_RescaleAndModify; 439 PetscCall(PetscMalloc1(10, &gl->schemes)); 440 gl->nschemes = 0; 441 442 { 443 /* p=1,q=1, r=s=2, A- and L-stable with error estimates of order 2 and 3 444 * Listed in Butcher & Podhaisky 2006. On error estimation in general linear methods for stiff ODE. 445 * irks(0.3,0,[.3,1],[1],1) 446 * Note: can be made to have classical order (not stage order) 2 by replacing 0.3 with 1-sqrt(1/2) 447 * but doing so would sacrifice the error estimator. 448 */ 449 const PetscScalar c[2] = {3. / 10., 1.}; 450 const PetscScalar a[2][2] = { 451 {3. / 10., 0 }, 452 {7. / 10., 3. / 10.} 453 }; 454 const PetscScalar b[2][2] = { 455 {7. / 10., 3. / 10.}, 456 {0, 1 } 457 }; 458 const PetscScalar u[2][2] = { 459 {1, 0}, 460 {1, 0} 461 }; 462 const PetscScalar v[2][2] = { 463 {1, 0}, 464 {0, 0} 465 }; 466 PetscCall(TSGLLESchemeCreate(1, 1, 2, 2, c, *a, *b, *u, *v, &gl->schemes[gl->nschemes++])); 467 } 468 469 { 470 /* p=q=2, r=s=3: irks(4/9,0,[1:3]/3,[0.33852],1) */ 471 /* http://www.math.auckland.ac.nz/~hpod/atlas/i2a.html */ 472 const PetscScalar c[3] = {1. / 3., 2. / 3., 1}; 473 const PetscScalar a[3][3] = { 474 {4. / 9., 0, 0 }, 475 {1.03750643704090e+00, 4. / 9., 0 }, 476 {7.67024779410304e-01, -3.81140216918943e-01, 4. / 9.} 477 }; 478 const PetscScalar b[3][3] = { 479 {0.767024779410304, -0.381140216918943, 4. / 9. }, 480 {0.000000000000000, 0.000000000000000, 1.000000000000000}, 481 {-2.075048385225385, 0.621728385225383, 1.277197204924873} 482 }; 483 const PetscScalar u[3][3] = { 484 {1.0000000000000000, -0.1111111111111109, -0.0925925925925922}, 485 {1.0000000000000000, -0.8152842148186744, -0.4199095530877056}, 486 {1.0000000000000000, 0.1696709930641948, 0.0539741070314165 } 487 }; 488 const PetscScalar v[3][3] = { 489 {1.0000000000000000, 0.1696709930641948, 0.0539741070314165}, 490 {0.000000000000000, 0.000000000000000, 0.000000000000000 }, 491 {0.000000000000000, 0.176122795075129, 0.000000000000000 } 492 }; 493 PetscCall(TSGLLESchemeCreate(2, 2, 3, 3, c, *a, *b, *u, *v, &gl->schemes[gl->nschemes++])); 494 } 495 { 496 /* p=q=3, r=s=4: irks(9/40,0,[1:4]/4,[0.3312 1.0050],[0.49541 1;1 0]) */ 497 const PetscScalar c[4] = {0.25, 0.5, 0.75, 1.0}; 498 const PetscScalar a[4][4] = { 499 {9. / 40., 0, 0, 0 }, 500 {2.11286958887701e-01, 9. / 40., 0, 0 }, 501 {9.46338294287584e-01, -3.42942861246094e-01, 9. / 40., 0 }, 502 {0.521490453970721, -0.662474225622980, 0.490476425459734, 9. / 40.} 503 }; 504 const PetscScalar b[4][4] = { 505 {0.521490453970721, -0.662474225622980, 0.490476425459734, 9. / 40. }, 506 {0.000000000000000, 0.000000000000000, 0.000000000000000, 1.000000000000000}, 507 {-0.084677029310348, 1.390757514776085, -1.568157386206001, 2.023192696767826}, 508 {0.465383797936408, 1.478273530625148, -1.930836081010182, 1.644872111193354} 509 }; 510 const PetscScalar u[4][4] = { 511 {1.00000000000000000, 0.02500000000001035, -0.02499999999999053, -0.00442708333332865}, 512 {1.00000000000000000, 0.06371304111232945, -0.04032173972189845, -0.01389438413189452}, 513 {1.00000000000000000, -0.07839543304147778, 0.04738685705116663, 0.02032603595928376 }, 514 {1.00000000000000000, 0.42550734619251651, 0.10800718022400080, -0.01726712647760034} 515 }; 516 const PetscScalar v[4][4] = { 517 {1.00000000000000000, 0.42550734619251651, 0.10800718022400080, -0.01726712647760034}, 518 {0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000 }, 519 {0.000000000000000, -1.761115796027561, -0.521284157173780, 0.258249384305463 }, 520 {0.000000000000000, -1.657693358744728, -1.052227765232394, 0.521284157173780 } 521 }; 522 PetscCall(TSGLLESchemeCreate(3, 3, 4, 4, c, *a, *b, *u, *v, &gl->schemes[gl->nschemes++])); 523 } 524 { 525 /* p=q=4, r=s=5: 526 irks(3/11,0,[1:5]/5, [0.1715 -0.1238 0.6617],... 527 [ -0.0812 0.4079 1.0000 528 1.0000 0 0 529 0.8270 1.0000 0]) 530 */ 531 const PetscScalar c[5] = {0.2, 0.4, 0.6, 0.8, 1.0}; 532 const PetscScalar a[5][5] = { 533 {2.72727272727352e-01, 0.00000000000000e+00, 0.00000000000000e+00, 0.00000000000000e+00, 0.00000000000000e+00}, 534 {-1.03980153733431e-01, 2.72727272727405e-01, 0.00000000000000e+00, 0.00000000000000e+00, 0.00000000000000e+00}, 535 {-1.58615400341492e+00, 7.44168951881122e-01, 2.72727272727309e-01, 0.00000000000000e+00, 0.00000000000000e+00}, 536 {-8.73658042865628e-01, 5.37884671894595e-01, -1.63298538799523e-01, 2.72727272726996e-01, 0.00000000000000e+00}, 537 {2.95489397443992e-01, -1.18481693910097e+00, -6.68029812659953e-01, 1.00716687860943e+00, 2.72727272727288e-01} 538 }; 539 const PetscScalar b[5][5] = { 540 {2.95489397443992e-01, -1.18481693910097e+00, -6.68029812659953e-01, 1.00716687860943e+00, 2.72727272727288e-01}, 541 {0.00000000000000e+00, 1.11022302462516e-16, -2.22044604925031e-16, 0.00000000000000e+00, 1.00000000000000e+00}, 542 {-4.05882503986005e+00, -4.00924006567769e+00, -1.38930610972481e+00, 4.45223930308488e+00, 6.32331093108427e-01}, 543 {8.35690179937017e+00, -2.26640927349732e+00, 6.86647884973826e+00, -5.22595158025740e+00, 4.50893068837431e+00}, 544 {1.27656267027479e+01, 2.80882153840821e+00, 8.91173096522890e+00, -1.07936444078906e+01, 4.82534148988854e+00} 545 }; 546 const PetscScalar u[5][5] = { 547 {1.00000000000000e+00, -7.27272727273551e-02, -3.45454545454419e-02, -4.12121212119565e-03, -2.96969696964014e-04}, 548 {1.00000000000000e+00, 2.31252881006154e-01, -8.29487834416481e-03, -9.07191207681020e-03, -1.70378403743473e-03}, 549 {1.00000000000000e+00, 1.16925777880663e+00, 3.59268562942635e-02, -4.09013451730615e-02, -1.02411119670164e-02}, 550 {1.00000000000000e+00, 1.02634463704356e+00, 1.59375044913405e-01, 1.89673015035370e-03, -4.89987231897569e-03}, 551 {1.00000000000000e+00, 1.27746320298021e+00, 2.37186008132728e-01, -8.28694373940065e-02, -5.34396510196430e-02} 552 }; 553 const PetscScalar v[5][5] = { 554 {1.00000000000000e+00, 1.27746320298021e+00, 2.37186008132728e-01, -8.28694373940065e-02, -5.34396510196430e-02}, 555 {0.00000000000000e+00, -1.77635683940025e-15, -1.99840144432528e-15, -9.99200722162641e-16, -3.33066907387547e-16}, 556 {0.00000000000000e+00, 4.37280081906924e+00, 5.49221645016377e-02, -8.88913177394943e-02, 1.12879077989154e-01 }, 557 {0.00000000000000e+00, -1.22399504837280e+01, -5.21287338448645e+00, -8.03952325565291e-01, 4.60298678047147e-01 }, 558 {0.00000000000000e+00, -1.85178762883829e+01, -5.21411849862624e+00, -1.04283436528809e+00, 7.49030161063651e-01 } 559 }; 560 PetscCall(TSGLLESchemeCreate(4, 4, 5, 5, c, *a, *b, *u, *v, &gl->schemes[gl->nschemes++])); 561 } 562 { 563 /* p=q=5, r=s=6; 564 irks(1/3,0,[1:6]/6,... 565 [-0.0489 0.4228 -0.8814 0.9021],... 566 [-0.3474 -0.6617 0.6294 0.2129 567 0.0044 -0.4256 -0.1427 -0.8936 568 -0.8267 0.4821 0.1371 -0.2557 569 -0.4426 -0.3855 -0.7514 0.3014]) 570 */ 571 const PetscScalar c[6] = {1. / 6, 2. / 6, 3. / 6, 4. / 6, 5. / 6, 1.}; 572 const PetscScalar a[6][6] = { 573 {3.33333333333940e-01, 0, 0, 0, 0, 0 }, 574 {-8.64423857333350e-02, 3.33333333332888e-01, 0, 0, 0, 0 }, 575 {-2.16850174258252e+00, -2.23619072028839e+00, 3.33333333335204e-01, 0, 0, 0 }, 576 {-4.73160970138997e+00, -3.89265344629268e+00, -2.76318716520933e-01, 3.33333333335759e-01, 0, 0 }, 577 {-6.75187540297338e+00, -7.90756533769377e+00, 7.90245051802259e-01, -4.48352364517632e-01, 3.33333333328483e-01, 0 }, 578 {-4.26488287921548e+00, -1.19320395589302e+01, 3.38924509887755e+00, -2.23969848002481e+00, 6.62807710124007e-01, 3.33333333335440e-01} 579 }; 580 const PetscScalar b[6][6] = { 581 {-4.26488287921548e+00, -1.19320395589302e+01, 3.38924509887755e+00, -2.23969848002481e+00, 6.62807710124007e-01, 3.33333333335440e-01 }, 582 {-8.88178419700125e-16, 4.44089209850063e-16, -1.54737334057131e-15, -8.88178419700125e-16, 0.00000000000000e+00, 1.00000000000001e+00 }, 583 {-2.87780425770651e+01, -1.13520448264971e+01, 2.62002318943161e+01, 2.56943874812797e+01, -3.06702268304488e+01, 6.68067773510103e+00 }, 584 {5.47971245256474e+01, 6.80366875868284e+01, -6.50952588861999e+01, -8.28643975339097e+01, 8.17416943896414e+01, -1.17819043489036e+01}, 585 {-2.33332114788869e+02, 6.12942539462634e+01, -4.91850135865944e+01, 1.82716844135480e+02, -1.29788173979395e+02, 3.09968095651099e+01 }, 586 {-1.72049132343751e+02, 8.60194713593999e+00, 7.98154219170200e-01, 1.50371386053218e+02, -1.18515423962066e+02, 2.50898277784663e+01 } 587 }; 588 const PetscScalar u[6][6] = { 589 {1.00000000000000e+00, -1.66666666666870e-01, -4.16666666664335e-02, -3.85802469124815e-03, -2.25051440302250e-04, -9.64506172339142e-06}, 590 {1.00000000000000e+00, 8.64423857327162e-02, -4.11484912671353e-02, -1.11450903217645e-02, -1.47651050487126e-03, -1.34395070766826e-04}, 591 {1.00000000000000e+00, 4.57135912953434e+00, 1.06514719719137e+00, 1.33517564218007e-01, 1.11365952968659e-02, 6.12382756769504e-04 }, 592 {1.00000000000000e+00, 9.23391519753404e+00, 2.22431212392095e+00, 2.91823807741891e-01, 2.52058456411084e-02, 1.22800542949647e-03 }, 593 {1.00000000000000e+00, 1.48175480533865e+01, 3.73439117461835e+00, 5.14648336541804e-01, 4.76430038853402e-02, 2.56798515502156e-03 }, 594 {1.00000000000000e+00, 1.50512347758335e+01, 4.10099701165164e+00, 5.66039141003603e-01, 3.91213893800891e-02, -2.99136269067853e-03} 595 }; 596 const PetscScalar v[6][6] = { 597 {1.00000000000000e+00, 1.50512347758335e+01, 4.10099701165164e+00, 5.66039141003603e-01, 3.91213893800891e-02, -2.99136269067853e-03}, 598 {0.00000000000000e+00, -4.88498130835069e-15, -6.43929354282591e-15, -3.55271367880050e-15, -1.22124532708767e-15, -3.12250225675825e-16}, 599 {0.00000000000000e+00, 1.22250171233141e+01, -1.77150760606169e+00, 3.54516769879390e-01, 6.22298845883398e-01, 2.31647447450276e-01 }, 600 {0.00000000000000e+00, -4.48339457331040e+01, -3.57363126641880e-01, 5.18750173123425e-01, 6.55727990241799e-02, 1.63175368287079e-01 }, 601 {0.00000000000000e+00, 1.37297394708005e+02, -1.60145272991317e+00, -5.05319555199441e+00, 1.55328940390990e-01, 9.16629423682464e-01 }, 602 {0.00000000000000e+00, 1.05703241119022e+02, -1.16610260983038e+00, -2.99767252773859e+00, -1.13472315553890e-01, 1.09742849254729e+00 } 603 }; 604 PetscCall(TSGLLESchemeCreate(5, 5, 6, 6, c, *a, *b, *u, *v, &gl->schemes[gl->nschemes++])); 605 } 606 PetscFunctionReturn(0); 607 } 608 609 /*@C 610 TSGLLESetType - sets the class of general linear method to use for time-stepping 611 612 Collective on TS 613 614 Input Parameters: 615 + ts - the TS context 616 - type - a method 617 618 Options Database Key: 619 . -ts_gl_type <type> - sets the method, use -help for a list of available method (e.g. irks) 620 621 Notes: 622 See "petsc/include/petscts.h" for available methods (for instance) 623 . TSGLLE_IRKS - Diagonally implicit methods with inherent Runge-Kutta stability (for stiff problems) 624 625 Normally, it is best to use the TSSetFromOptions() command and 626 then set the TSGLLE type from the options database rather than by using 627 this routine. Using the options database provides the user with 628 maximum flexibility in evaluating the many different solvers. 629 The TSGLLESetType() routine is provided for those situations where it 630 is necessary to set the timestepping solver independently of the 631 command line or options database. This might be the case, for example, 632 when the choice of solver changes during the execution of the 633 program, and the user's application is taking responsibility for 634 choosing the appropriate method. 635 636 Level: intermediate 637 638 @*/ 639 PetscErrorCode TSGLLESetType(TS ts, TSGLLEType type) 640 { 641 PetscFunctionBegin; 642 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 643 PetscValidCharPointer(type, 2); 644 PetscTryMethod(ts, "TSGLLESetType_C", (TS, TSGLLEType), (ts, type)); 645 PetscFunctionReturn(0); 646 } 647 648 /*@C 649 TSGLLESetAcceptType - sets the acceptance test 650 651 Time integrators that need to control error must have the option to reject a time step based on local error 652 estimates. This function allows different schemes to be set. 653 654 Logically Collective on TS 655 656 Input Parameters: 657 + ts - the TS context 658 - type - the type 659 660 Options Database Key: 661 . -ts_gl_accept_type <type> - sets the method used to determine whether to accept or reject a step 662 663 Level: intermediate 664 665 .seealso: `TS`, `TSGLLE`, `TSGLLEAcceptRegister()`, `TSGLLEAdapt`, `set` `type` 666 @*/ 667 PetscErrorCode TSGLLESetAcceptType(TS ts, TSGLLEAcceptType type) 668 { 669 PetscFunctionBegin; 670 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 671 PetscValidCharPointer(type, 2); 672 PetscTryMethod(ts, "TSGLLESetAcceptType_C", (TS, TSGLLEAcceptType), (ts, type)); 673 PetscFunctionReturn(0); 674 } 675 676 /*@C 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 Notes: 688 This allows the user set options on the TSGLLEAdapt object. Usually it is better to do this using the options 689 database, so this function is rarely needed. 690 691 Level: advanced 692 693 .seealso: `TSGLLEAdapt`, `TSGLLEAdaptRegister()` 694 @*/ 695 PetscErrorCode TSGLLEGetAdapt(TS ts, TSGLLEAdapt *adapt) 696 { 697 PetscFunctionBegin; 698 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 699 PetscValidPointer(adapt, 2); 700 PetscUseMethod(ts, "TSGLLEGetAdapt_C", (TS, TSGLLEAdapt *), (ts, adapt)); 701 PetscFunctionReturn(0); 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(0); 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(0); 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 PetscCall(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(0); 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(0); 757 PetscCall((*gl->Destroy)(gl)); 758 } 759 760 PetscCall(PetscFunctionListFind(TSGLLEList, type, &r)); 761 PetscCheck(r, PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown TSGLLE type \"%s\" given", type); 762 PetscCall((*r)(ts)); 763 PetscCall(PetscStrcpy(gl->type_name, type)); 764 PetscFunctionReturn(0); 765 } 766 767 static PetscErrorCode TSGLLESetAcceptType_GLLE(TS ts, TSGLLEAcceptType type) 768 { 769 TSGLLEAcceptFunction r; 770 TS_GLLE *gl = (TS_GLLE *)ts->data; 771 772 PetscFunctionBegin; 773 PetscCall(PetscFunctionListFind(TSGLLEAcceptList, type, &r)); 774 PetscCheck(r, PETSC_COMM_SELF, 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(0); 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(0); 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(0); 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(0); 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 > 0 && ++ts->num_snes_failures >= ts->max_snes_failures) { 865 ts->reason = TS_DIVERGED_NONLINEAR_SOLVE; 866 PetscCall(PetscInfo(ts, "Step=%" PetscInt_FMT ", nonlinear solve solve failures %" PetscInt_FMT " greater than current TS allowed, stopping solve\n", ts->steps, ts->num_snes_failures)); 867 PetscFunctionReturn(0); 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 > 0 && ++ts->num_snes_failures >= ts->max_snes_failures) { 940 ts->reason = TS_DIVERGED_NONLINEAR_SOLVE; 941 PetscCall(PetscInfo(ts, "Step=%" PetscInt_FMT ", nonlinear solve solve failures %" PetscInt_FMT " greater than current TS allowed, stopping solve\n", ts->steps, ts->num_snes_failures)); 942 PetscFunctionReturn(0); 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(0); 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(0); 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(0); 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(0); 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(0); 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(0); 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(PETSC_COMM_SELF, 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(0); 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(0); 1195 } 1196 1197 /*@C 1198 TSGLLERegister - adds a TSGLLE implementation 1199 1200 Not Collective 1201 1202 Input Parameters: 1203 + name_scheme - name of user-defined general linear scheme 1204 - routine_create - routine to create method context 1205 1206 Notes: 1207 TSGLLERegister() may be called multiple times to add several user-defined families. 1208 1209 Sample usage: 1210 .vb 1211 TSGLLERegister("my_scheme",MySchemeCreate); 1212 .ve 1213 1214 Then, your scheme can be chosen with the procedural interface via 1215 $ TSGLLESetType(ts,"my_scheme") 1216 or at runtime via the option 1217 $ -ts_gl_type my_scheme 1218 1219 Level: advanced 1220 1221 .seealso: `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(0); 1229 } 1230 1231 /*@C 1232 TSGLLEAcceptRegister - adds a TSGLLE acceptance scheme 1233 1234 Not Collective 1235 1236 Input Parameters: 1237 + name_scheme - name of user-defined acceptance scheme 1238 - routine_create - routine to create method context 1239 1240 Notes: 1241 TSGLLEAcceptRegister() may be called multiple times to add several user-defined families. 1242 1243 Sample usage: 1244 .vb 1245 TSGLLEAcceptRegister("my_scheme",MySchemeCreate); 1246 .ve 1247 1248 Then, your scheme can be chosen with the procedural interface via 1249 $ TSGLLESetAcceptType(ts,"my_scheme") 1250 or at runtime via the option 1251 $ -ts_gl_accept_type my_scheme 1252 1253 Level: advanced 1254 1255 .seealso: `TSGLLERegisterAll()` 1256 @*/ 1257 PetscErrorCode TSGLLEAcceptRegister(const char sname[], TSGLLEAcceptFunction function) 1258 { 1259 PetscFunctionBegin; 1260 PetscCall(PetscFunctionListAdd(&TSGLLEAcceptList, sname, function)); 1261 PetscFunctionReturn(0); 1262 } 1263 1264 /*@C 1265 TSGLLERegisterAll - Registers all of the general linear methods in TSGLLE 1266 1267 Not Collective 1268 1269 Level: advanced 1270 1271 .seealso: `TSGLLERegisterDestroy()` 1272 @*/ 1273 PetscErrorCode TSGLLERegisterAll(void) 1274 { 1275 PetscFunctionBegin; 1276 if (TSGLLERegisterAllCalled) PetscFunctionReturn(0); 1277 TSGLLERegisterAllCalled = PETSC_TRUE; 1278 1279 PetscCall(TSGLLERegister(TSGLLE_IRKS, TSGLLECreate_IRKS)); 1280 PetscCall(TSGLLEAcceptRegister(TSGLLEACCEPT_ALWAYS, TSGLLEAccept_Always)); 1281 PetscFunctionReturn(0); 1282 } 1283 1284 /*@C 1285 TSGLLEInitializePackage - This function initializes everything in the TSGLLE package. It is called 1286 from TSInitializePackage(). 1287 1288 Level: developer 1289 1290 .seealso: `PetscInitialize()` 1291 @*/ 1292 PetscErrorCode TSGLLEInitializePackage(void) 1293 { 1294 PetscFunctionBegin; 1295 if (TSGLLEPackageInitialized) PetscFunctionReturn(0); 1296 TSGLLEPackageInitialized = PETSC_TRUE; 1297 PetscCall(TSGLLERegisterAll()); 1298 PetscCall(PetscRegisterFinalize(TSGLLEFinalizePackage)); 1299 PetscFunctionReturn(0); 1300 } 1301 1302 /*@C 1303 TSGLLEFinalizePackage - This function destroys everything in the TSGLLE package. It is 1304 called from PetscFinalize(). 1305 1306 Level: developer 1307 1308 .seealso: `PetscFinalize()` 1309 @*/ 1310 PetscErrorCode TSGLLEFinalizePackage(void) 1311 { 1312 PetscFunctionBegin; 1313 PetscCall(PetscFunctionListDestroy(&TSGLLEList)); 1314 PetscCall(PetscFunctionListDestroy(&TSGLLEAcceptList)); 1315 TSGLLEPackageInitialized = PETSC_FALSE; 1316 TSGLLERegisterAllCalled = PETSC_FALSE; 1317 PetscFunctionReturn(0); 1318 } 1319 1320 /* ------------------------------------------------------------ */ 1321 /*MC 1322 TSGLLE - DAE solver using implicit General Linear methods 1323 1324 These methods contain Runge-Kutta and multistep schemes as special cases. These special cases have some fundamental 1325 limitations. For example, diagonally implicit Runge-Kutta cannot have stage order greater than 1 which limits their 1326 applicability to very stiff systems. Meanwhile, multistep methods cannot be A-stable for order greater than 2 and BDF 1327 are not 0-stable for order greater than 6. GL methods can be A- and L-stable with arbitrarily high stage order and 1328 reliable error estimates for both 1 and 2 orders higher to facilitate adaptive step sizes and adaptive order schemes. 1329 All this is possible while preserving a singly diagonally implicit structure. 1330 1331 Options database keys: 1332 + -ts_gl_type <type> - the class of general linear method (irks) 1333 . -ts_gl_rtol <tol> - relative error 1334 . -ts_gl_atol <tol> - absolute error 1335 . -ts_gl_min_order <p> - minimum order method to consider (default=1) 1336 . -ts_gl_max_order <p> - maximum order method to consider (default=3) 1337 . -ts_gl_start_order <p> - order of starting method (default=1) 1338 . -ts_gl_complete <method> - method to use for completing the step (rescale-and-modify or rescale) 1339 - -ts_adapt_type <method> - adaptive controller to use (none step both) 1340 1341 Notes: 1342 This integrator can be applied to DAE. 1343 1344 Diagonally implicit general linear (DIGL) methods are a generalization of diagonally implicit Runge-Kutta (DIRK). 1345 They are represented by the tableau 1346 1347 .vb 1348 A | U 1349 ------- 1350 B | V 1351 .ve 1352 1353 combined with a vector c of abscissa. "Diagonally implicit" means that A is lower triangular. 1354 A step of the general method reads 1355 1356 .vb 1357 [ Y ] = [A U] [ Y' ] 1358 [X^k] = [B V] [X^{k-1}] 1359 .ve 1360 1361 where Y is the multivector of stage values, Y' is the multivector of stage derivatives, X^k is the Nordsieck vector of 1362 the solution at step k. The Nordsieck vector consists of the first r moments of the solution, given by 1363 1364 .vb 1365 X = [x_0,x_1,...,x_{r-1}] = [x, h x', h^2 x'', ..., h^{r-1} x^{(r-1)} ] 1366 .ve 1367 1368 If A is lower triangular, we can solve the stages (Y,Y') sequentially 1369 1370 .vb 1371 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} 1372 .ve 1373 1374 and then construct the pieces to carry to the next step 1375 1376 .vb 1377 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} 1378 .ve 1379 1380 Note that when the equations are cast in implicit form, we are using the stage equation to define y'_i 1381 in terms of y_i and known stuff (y_j for j<i and x_j for all j). 1382 1383 Error estimation 1384 1385 At present, the most attractive GL methods for stiff problems are singly diagonally implicit schemes which posses 1386 Inherent Runge-Kutta Stability (IRKS). These methods have r=s, the number of items passed between steps is equal to 1387 the number of stages. The order and stage-order are one less than the number of stages. We use the error estimates 1388 in the 2007 paper which provide the following estimates 1389 1390 .vb 1391 h^{p+1} X^{(p+1)} = phi_0^T Y' + [0 psi_0^T] Xold 1392 h^{p+2} X^{(p+2)} = phi_1^T Y' + [0 psi_1^T] Xold 1393 h^{p+2} (dx'/dx) X^{(p+1)} = phi_2^T Y' + [0 psi_2^T] Xold 1394 .ve 1395 1396 These estimates are accurate to O(h^{p+3}). 1397 1398 Changing the step size 1399 1400 We use the generalized "rescale and modify" scheme, see equation (4.5) of the 2007 paper. 1401 1402 Level: beginner 1403 1404 References: 1405 + * - John Butcher and Z. Jackieweicz and W. Wright, On error propagation in general linear methods for 1406 ordinary differential equations, Journal of Complexity, Vol 23, 2007. 1407 - * - John Butcher, Numerical methods for ordinary differential equations, second edition, Wiley, 2009. 1408 1409 .seealso: `TSCreate()`, `TS`, `TSSetType()` 1410 1411 M*/ 1412 PETSC_EXTERN PetscErrorCode TSCreate_GLLE(TS ts) 1413 { 1414 TS_GLLE *gl; 1415 1416 PetscFunctionBegin; 1417 PetscCall(TSGLLEInitializePackage()); 1418 1419 PetscCall(PetscNew(&gl)); 1420 ts->data = (void *)gl; 1421 1422 ts->ops->reset = TSReset_GLLE; 1423 ts->ops->destroy = TSDestroy_GLLE; 1424 ts->ops->view = TSView_GLLE; 1425 ts->ops->setup = TSSetUp_GLLE; 1426 ts->ops->solve = TSSolve_GLLE; 1427 ts->ops->setfromoptions = TSSetFromOptions_GLLE; 1428 ts->ops->snesfunction = SNESTSFormFunction_GLLE; 1429 ts->ops->snesjacobian = SNESTSFormJacobian_GLLE; 1430 1431 ts->usessnes = PETSC_TRUE; 1432 1433 gl->max_step_rejections = 1; 1434 gl->min_order = 1; 1435 gl->max_order = 3; 1436 gl->start_order = 1; 1437 gl->current_scheme = -1; 1438 gl->extrapolate = PETSC_FALSE; 1439 1440 gl->wrms_atol = 1e-8; 1441 gl->wrms_rtol = 1e-5; 1442 1443 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSGLLESetType_C", &TSGLLESetType_GLLE)); 1444 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSGLLESetAcceptType_C", &TSGLLESetAcceptType_GLLE)); 1445 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSGLLEGetAdapt_C", &TSGLLEGetAdapt_GLLE)); 1446 PetscFunctionReturn(0); 1447 } 1448