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(PETSC_SUCCESS); 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(PETSC_SUCCESS); 61 } 62 63 static PetscErrorCode DMCoarsenHook_TSGLLE(DM fine, DM coarse, void *ctx) 64 { 65 PetscFunctionBegin; 66 PetscFunctionReturn(PETSC_SUCCESS); 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(PETSC_SUCCESS); 82 } 83 84 static PetscErrorCode DMSubDomainHook_TSGLLE(DM dm, DM subdm, void *ctx) 85 { 86 PetscFunctionBegin; 87 PetscFunctionReturn(PETSC_SUCCESS); 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(PETSC_SUCCESS); 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(PETSC_SUCCESS); 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(PETSC_SUCCESS); 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(PETSC_SUCCESS); 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(PETSC_SUCCESS); 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(PETSC_SUCCESS); 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(PETSC_SUCCESS); 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(PETSC_SUCCESS); 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(PETSC_SUCCESS); 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(PETSC_SUCCESS); 607 } 608 609 /*@C 610 TSGLLESetType - sets the class of general linear method, `TSGLLE` to use for time-stepping 611 612 Collective 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 Level: intermediate 622 623 Notes: 624 See "petsc/include/petscts.h" for available methods (for instance) 625 . TSGLLE_IRKS - Diagonally implicit methods with inherent Runge-Kutta stability (for stiff problems) 626 627 Normally, it is best to use the `TSSetFromOptions()` command and 628 then set the `TSGLLE` type from the options database rather than by using 629 this routine. Using the options database provides the user with 630 maximum flexibility in evaluating the many different solvers. 631 The `TSGLLESetType()` routine is provided for those situations where it 632 is necessary to set the timestepping solver independently of the 633 command line or options database. This might be the case, for example, 634 when the choice of solver changes during the execution of the 635 program, and the user's application is taking responsibility for 636 choosing the appropriate method. 637 638 .seealso: [](chapter_ts), `TS`, `TSGLLEType`, `TSGLLE` 639 @*/ 640 PetscErrorCode TSGLLESetType(TS ts, TSGLLEType type) 641 { 642 PetscFunctionBegin; 643 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 644 PetscValidCharPointer(type, 2); 645 PetscTryMethod(ts, "TSGLLESetType_C", (TS, TSGLLEType), (ts, type)); 646 PetscFunctionReturn(PETSC_SUCCESS); 647 } 648 649 /*@C 650 TSGLLESetAcceptType - sets the acceptance test for `TSGLLE` 651 652 Time integrators that need to control error must have the option to reject a time step based on local error 653 estimates. This function allows different schemes to be set. 654 655 Logically Collective 656 657 Input Parameters: 658 + ts - the `TS` context 659 - type - the type 660 661 Options Database Key: 662 . -ts_gl_accept_type <type> - sets the method used to determine whether to accept or reject a step 663 664 Level: intermediate 665 666 .seealso: [](chapter_ts), `TS`, `TSGLLE`, `TSGLLEAcceptRegister()`, `TSGLLEAdapt` 667 @*/ 668 PetscErrorCode TSGLLESetAcceptType(TS ts, TSGLLEAcceptType type) 669 { 670 PetscFunctionBegin; 671 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 672 PetscValidCharPointer(type, 2); 673 PetscTryMethod(ts, "TSGLLESetAcceptType_C", (TS, TSGLLEAcceptType), (ts, type)); 674 PetscFunctionReturn(PETSC_SUCCESS); 675 } 676 677 /*@C 678 TSGLLEGetAdapt - gets the `TSGLLEAdapt` object from the `TS` 679 680 Not Collective 681 682 Input Parameter: 683 . ts - the `TS` context 684 685 Output Parameter: 686 . adapt - the `TSGLLEAdapt` context 687 688 Level: advanced 689 690 Note: 691 This allows the user set options on the `TSGLLEAdapt` object. Usually it is better to do this using the options 692 database, so this function is rarely needed. 693 694 .seealso: [](chapter_ts), `TS`, `TSGLLE`, `TSGLLEAdapt`, `TSGLLEAdaptRegister()` 695 @*/ 696 PetscErrorCode TSGLLEGetAdapt(TS ts, TSGLLEAdapt *adapt) 697 { 698 PetscFunctionBegin; 699 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 700 PetscValidPointer(adapt, 2); 701 PetscUseMethod(ts, "TSGLLEGetAdapt_C", (TS, TSGLLEAdapt *), (ts, adapt)); 702 PetscFunctionReturn(PETSC_SUCCESS); 703 } 704 705 static PetscErrorCode TSGLLEAccept_Always(TS ts, PetscReal tleft, PetscReal h, const PetscReal enorms[], PetscBool *accept) 706 { 707 PetscFunctionBegin; 708 *accept = PETSC_TRUE; 709 PetscFunctionReturn(PETSC_SUCCESS); 710 } 711 712 static PetscErrorCode TSGLLEUpdateWRMS(TS ts) 713 { 714 TS_GLLE *gl = (TS_GLLE *)ts->data; 715 PetscScalar *x, *w; 716 PetscInt n, i; 717 718 PetscFunctionBegin; 719 PetscCall(VecGetArray(gl->X[0], &x)); 720 PetscCall(VecGetArray(gl->W, &w)); 721 PetscCall(VecGetLocalSize(gl->W, &n)); 722 for (i = 0; i < n; i++) w[i] = 1. / (gl->wrms_atol + gl->wrms_rtol * PetscAbsScalar(x[i])); 723 PetscCall(VecRestoreArray(gl->X[0], &x)); 724 PetscCall(VecRestoreArray(gl->W, &w)); 725 PetscFunctionReturn(PETSC_SUCCESS); 726 } 727 728 static PetscErrorCode TSGLLEVecNormWRMS(TS ts, Vec X, PetscReal *nrm) 729 { 730 TS_GLLE *gl = (TS_GLLE *)ts->data; 731 PetscScalar *x, *w; 732 PetscReal sum = 0.0, gsum; 733 PetscInt n, N, i; 734 735 PetscFunctionBegin; 736 PetscCall(VecGetArray(X, &x)); 737 PetscCall(VecGetArray(gl->W, &w)); 738 PetscCall(VecGetLocalSize(gl->W, &n)); 739 for (i = 0; i < n; i++) sum += PetscAbsScalar(PetscSqr(x[i] * w[i])); 740 PetscCall(VecRestoreArray(X, &x)); 741 PetscCall(VecRestoreArray(gl->W, &w)); 742 PetscCall(MPIU_Allreduce(&sum, &gsum, 1, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)ts))); 743 PetscCall(VecGetSize(gl->W, &N)); 744 *nrm = PetscSqrtReal(gsum / (1. * N)); 745 PetscFunctionReturn(PETSC_SUCCESS); 746 } 747 748 static PetscErrorCode TSGLLESetType_GLLE(TS ts, TSGLLEType type) 749 { 750 PetscBool same; 751 TS_GLLE *gl = (TS_GLLE *)ts->data; 752 PetscErrorCode (*r)(TS); 753 754 PetscFunctionBegin; 755 if (gl->type_name[0]) { 756 PetscCall(PetscStrcmp(gl->type_name, type, &same)); 757 if (same) PetscFunctionReturn(PETSC_SUCCESS); 758 PetscCall((*gl->Destroy)(gl)); 759 } 760 761 PetscCall(PetscFunctionListFind(TSGLLEList, type, &r)); 762 PetscCheck(r, PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown TSGLLE type \"%s\" given", type); 763 PetscCall((*r)(ts)); 764 PetscCall(PetscStrncpy(gl->type_name, type, sizeof(gl->type_name))); 765 PetscFunctionReturn(PETSC_SUCCESS); 766 } 767 768 static PetscErrorCode TSGLLESetAcceptType_GLLE(TS ts, TSGLLEAcceptType type) 769 { 770 TSGLLEAcceptFunction r; 771 TS_GLLE *gl = (TS_GLLE *)ts->data; 772 773 PetscFunctionBegin; 774 PetscCall(PetscFunctionListFind(TSGLLEAcceptList, type, &r)); 775 PetscCheck(r, PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown TSGLLEAccept type \"%s\" given", type); 776 gl->Accept = r; 777 PetscCall(PetscStrncpy(gl->accept_name, type, sizeof(gl->accept_name))); 778 PetscFunctionReturn(PETSC_SUCCESS); 779 } 780 781 static PetscErrorCode TSGLLEGetAdapt_GLLE(TS ts, TSGLLEAdapt *adapt) 782 { 783 TS_GLLE *gl = (TS_GLLE *)ts->data; 784 785 PetscFunctionBegin; 786 if (!gl->adapt) { 787 PetscCall(TSGLLEAdaptCreate(PetscObjectComm((PetscObject)ts), &gl->adapt)); 788 PetscCall(PetscObjectIncrementTabLevel((PetscObject)gl->adapt, (PetscObject)ts, 1)); 789 } 790 *adapt = gl->adapt; 791 PetscFunctionReturn(PETSC_SUCCESS); 792 } 793 794 static PetscErrorCode TSGLLEChooseNextScheme(TS ts, PetscReal h, const PetscReal hmnorm[], PetscInt *next_scheme, PetscReal *next_h, PetscBool *finish) 795 { 796 TS_GLLE *gl = (TS_GLLE *)ts->data; 797 PetscInt i, n, cur_p, cur, next_sc, candidates[64], orders[64]; 798 PetscReal errors[64], costs[64], tleft; 799 800 PetscFunctionBegin; 801 cur = -1; 802 cur_p = gl->schemes[gl->current_scheme]->p; 803 tleft = ts->max_time - (ts->ptime + ts->time_step); 804 for (i = 0, n = 0; i < gl->nschemes; i++) { 805 TSGLLEScheme sc = gl->schemes[i]; 806 if (sc->p < gl->min_order || gl->max_order < sc->p) continue; 807 if (sc->p == cur_p - 1) errors[n] = PetscAbsScalar(sc->alpha[0]) * hmnorm[0]; 808 else if (sc->p == cur_p) errors[n] = PetscAbsScalar(sc->alpha[0]) * hmnorm[1]; 809 else if (sc->p == cur_p + 1) errors[n] = PetscAbsScalar(sc->alpha[0]) * (hmnorm[2] + hmnorm[3]); 810 else continue; 811 candidates[n] = i; 812 orders[n] = PetscMin(sc->p, sc->q); /* order of global truncation error */ 813 costs[n] = sc->s; /* estimate the cost as the number of stages */ 814 if (i == gl->current_scheme) cur = n; 815 n++; 816 } 817 PetscCheck(cur >= 0 && gl->nschemes > cur, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Current scheme not found in scheme list"); 818 PetscCall(TSGLLEAdaptChoose(gl->adapt, n, orders, errors, costs, cur, h, tleft, &next_sc, next_h, finish)); 819 *next_scheme = candidates[next_sc]; 820 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, 821 gl->schemes[*next_scheme]->r, gl->schemes[*next_scheme]->s, (double)*next_h, PetscBools[*finish])); 822 PetscFunctionReturn(PETSC_SUCCESS); 823 } 824 825 static PetscErrorCode TSGLLEGetMaxSizes(TS ts, PetscInt *max_r, PetscInt *max_s) 826 { 827 TS_GLLE *gl = (TS_GLLE *)ts->data; 828 829 PetscFunctionBegin; 830 *max_r = gl->schemes[gl->nschemes - 1]->r; 831 *max_s = gl->schemes[gl->nschemes - 1]->s; 832 PetscFunctionReturn(PETSC_SUCCESS); 833 } 834 835 static PetscErrorCode TSSolve_GLLE(TS ts) 836 { 837 TS_GLLE *gl = (TS_GLLE *)ts->data; 838 PetscInt i, k, its, lits, max_r, max_s; 839 PetscBool final_step, finish; 840 SNESConvergedReason snesreason; 841 842 PetscFunctionBegin; 843 PetscCall(TSMonitor(ts, ts->steps, ts->ptime, ts->vec_sol)); 844 845 PetscCall(TSGLLEGetMaxSizes(ts, &max_r, &max_s)); 846 PetscCall(VecCopy(ts->vec_sol, gl->X[0])); 847 for (i = 1; i < max_r; i++) PetscCall(VecZeroEntries(gl->X[i])); 848 PetscCall(TSGLLEUpdateWRMS(ts)); 849 850 if (0) { 851 /* Find consistent initial data for DAE */ 852 gl->stage_time = ts->ptime + ts->time_step; 853 gl->scoeff = 1.; 854 gl->stage = 0; 855 856 PetscCall(VecCopy(ts->vec_sol, gl->Z)); 857 PetscCall(VecCopy(ts->vec_sol, gl->Y)); 858 PetscCall(SNESSolve(ts->snes, NULL, gl->Y)); 859 PetscCall(SNESGetIterationNumber(ts->snes, &its)); 860 PetscCall(SNESGetLinearSolveIterations(ts->snes, &lits)); 861 PetscCall(SNESGetConvergedReason(ts->snes, &snesreason)); 862 863 ts->snes_its += its; 864 ts->ksp_its += lits; 865 if (snesreason < 0 && ts->max_snes_failures > 0 && ++ts->num_snes_failures >= ts->max_snes_failures) { 866 ts->reason = TS_DIVERGED_NONLINEAR_SOLVE; 867 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)); 868 PetscFunctionReturn(PETSC_SUCCESS); 869 } 870 } 871 872 PetscCheck(gl->current_scheme >= 0, PETSC_COMM_SELF, PETSC_ERR_ORDER, "A starting scheme has not been provided"); 873 874 for (k = 0, final_step = PETSC_FALSE, finish = PETSC_FALSE; k < ts->max_steps && !finish; k++) { 875 PetscInt j, r, s, next_scheme = 0, rejections; 876 PetscReal h, hmnorm[4], enorm[3], next_h; 877 PetscBool accept; 878 const PetscScalar *c, *a, *u; 879 Vec *X, *Ydot, Y; 880 TSGLLEScheme scheme = gl->schemes[gl->current_scheme]; 881 882 r = scheme->r; 883 s = scheme->s; 884 c = scheme->c; 885 a = scheme->a; 886 u = scheme->u; 887 h = ts->time_step; 888 X = gl->X; 889 Ydot = gl->Ydot; 890 Y = gl->Y; 891 892 if (ts->ptime > ts->max_time) break; 893 894 /* 895 We only call PreStep at the start of each STEP, not each STAGE. This is because it is 896 possible to fail (have to restart a step) after multiple stages. 897 */ 898 PetscCall(TSPreStep(ts)); 899 900 rejections = 0; 901 while (1) { 902 for (i = 0; i < s; i++) { 903 PetscScalar shift; 904 gl->scoeff = 1. / PetscRealPart(a[i * s + i]); 905 shift = gl->scoeff / ts->time_step; 906 gl->stage = i; 907 gl->stage_time = ts->ptime + PetscRealPart(c[i]) * h; 908 909 /* 910 * Stage equation: Y = h A Y' + U X 911 * We assume that A is lower-triangular so that we can solve the stages (Y,Y') sequentially 912 * Build the affine vector z_i = -[1/(h a_ii)](h sum_j a_ij y'_j + sum_j u_ij x_j) 913 * Then y'_i = z + 1/(h a_ii) y_i 914 */ 915 PetscCall(VecZeroEntries(gl->Z)); 916 for (j = 0; j < r; j++) PetscCall(VecAXPY(gl->Z, -shift * u[i * r + j], X[j])); 917 for (j = 0; j < i; j++) PetscCall(VecAXPY(gl->Z, -shift * h * a[i * s + j], Ydot[j])); 918 /* Note: Z is used within function evaluation, Ydot = Z + shift*Y */ 919 920 /* Compute an estimate of Y to start Newton iteration */ 921 if (gl->extrapolate) { 922 if (i == 0) { 923 /* Linear extrapolation on the first stage */ 924 PetscCall(VecWAXPY(Y, c[i] * h, X[1], X[0])); 925 } else { 926 /* Linear extrapolation from the last stage */ 927 PetscCall(VecAXPY(Y, (c[i] - c[i - 1]) * h, Ydot[i - 1])); 928 } 929 } else if (i == 0) { /* Directly use solution from the last step, otherwise reuse the last stage (do nothing) */ 930 PetscCall(VecCopy(X[0], Y)); 931 } 932 933 /* Solve this stage (Ydot[i] is computed during function evaluation) */ 934 PetscCall(SNESSolve(ts->snes, NULL, Y)); 935 PetscCall(SNESGetIterationNumber(ts->snes, &its)); 936 PetscCall(SNESGetLinearSolveIterations(ts->snes, &lits)); 937 PetscCall(SNESGetConvergedReason(ts->snes, &snesreason)); 938 ts->snes_its += its; 939 ts->ksp_its += lits; 940 if (snesreason < 0 && ts->max_snes_failures > 0 && ++ts->num_snes_failures >= ts->max_snes_failures) { 941 ts->reason = TS_DIVERGED_NONLINEAR_SOLVE; 942 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)); 943 PetscFunctionReturn(PETSC_SUCCESS); 944 } 945 } 946 947 gl->stage_time = ts->ptime + ts->time_step; 948 949 PetscCall((*gl->EstimateHigherMoments)(scheme, h, Ydot, gl->X, gl->himom)); 950 /* hmnorm[i] = h^{p+i}x^{(p+i)} with i=0,1,2; hmnorm[3] = h^{p+2}(dx'/dx) x^{(p+1)} */ 951 for (i = 0; i < 3; i++) PetscCall(TSGLLEVecNormWRMS(ts, gl->himom[i], &hmnorm[i + 1])); 952 enorm[0] = PetscRealPart(scheme->alpha[0]) * hmnorm[1]; 953 enorm[1] = PetscRealPart(scheme->beta[0]) * hmnorm[2]; 954 enorm[2] = PetscRealPart(scheme->gamma[0]) * hmnorm[3]; 955 PetscCall((*gl->Accept)(ts, ts->max_time - gl->stage_time, h, enorm, &accept)); 956 if (accept) goto accepted; 957 rejections++; 958 PetscCall(PetscInfo(ts, "Step %" PetscInt_FMT " (t=%g) not accepted, rejections=%" PetscInt_FMT "\n", k, (double)gl->stage_time, rejections)); 959 if (rejections > gl->max_step_rejections) break; 960 /* 961 There are lots of reasons why a step might be rejected, including solvers not converging and other factors that 962 TSGLLEChooseNextScheme does not support. Additionally, the error estimates may be very screwed up, so I'm not 963 convinced that it's safe to just compute a new error estimate using the same interface as the current adaptor 964 (the adaptor interface probably has to change). Here we make an arbitrary and naive choice. This assumes that 965 steps were written in Nordsieck form. The "correct" method would be to re-complete the previous time step with 966 the correct "next" step size. It is unclear to me whether the present ad-hoc method of rescaling X is stable. 967 */ 968 h *= 0.5; 969 for (i = 1; i < scheme->r; i++) PetscCall(VecScale(X[i], PetscPowRealInt(0.5, i))); 970 } 971 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); 972 973 accepted: 974 /* This term is not error, but it *would* be the leading term for a lower order method */ 975 PetscCall(TSGLLEVecNormWRMS(ts, gl->X[scheme->r - 1], &hmnorm[0])); 976 /* Correct scaling so that these are equivalent to norms of the Nordsieck vectors */ 977 978 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])); 979 if (!final_step) { 980 PetscCall(TSGLLEChooseNextScheme(ts, h, hmnorm, &next_scheme, &next_h, &final_step)); 981 } else { 982 /* Dummy values to complete the current step in a consistent manner */ 983 next_scheme = gl->current_scheme; 984 next_h = h; 985 finish = PETSC_TRUE; 986 } 987 988 X = gl->Xold; 989 gl->Xold = gl->X; 990 gl->X = X; 991 PetscCall((*gl->CompleteStep)(scheme, h, gl->schemes[next_scheme], next_h, Ydot, gl->Xold, gl->X)); 992 993 PetscCall(TSGLLEUpdateWRMS(ts)); 994 995 /* Post the solution for the user, we could avoid this copy with a small bit of cleverness */ 996 PetscCall(VecCopy(gl->X[0], ts->vec_sol)); 997 ts->ptime += h; 998 ts->steps++; 999 1000 PetscCall(TSPostEvaluate(ts)); 1001 PetscCall(TSPostStep(ts)); 1002 PetscCall(TSMonitor(ts, ts->steps, ts->ptime, ts->vec_sol)); 1003 1004 gl->current_scheme = next_scheme; 1005 ts->time_step = next_h; 1006 } 1007 PetscFunctionReturn(PETSC_SUCCESS); 1008 } 1009 1010 /*------------------------------------------------------------*/ 1011 1012 static PetscErrorCode TSReset_GLLE(TS ts) 1013 { 1014 TS_GLLE *gl = (TS_GLLE *)ts->data; 1015 PetscInt max_r, max_s; 1016 1017 PetscFunctionBegin; 1018 if (gl->setupcalled) { 1019 PetscCall(TSGLLEGetMaxSizes(ts, &max_r, &max_s)); 1020 PetscCall(VecDestroyVecs(max_r, &gl->Xold)); 1021 PetscCall(VecDestroyVecs(max_r, &gl->X)); 1022 PetscCall(VecDestroyVecs(max_s, &gl->Ydot)); 1023 PetscCall(VecDestroyVecs(3, &gl->himom)); 1024 PetscCall(VecDestroy(&gl->W)); 1025 PetscCall(VecDestroy(&gl->Y)); 1026 PetscCall(VecDestroy(&gl->Z)); 1027 } 1028 gl->setupcalled = PETSC_FALSE; 1029 PetscFunctionReturn(PETSC_SUCCESS); 1030 } 1031 1032 static PetscErrorCode TSDestroy_GLLE(TS ts) 1033 { 1034 TS_GLLE *gl = (TS_GLLE *)ts->data; 1035 1036 PetscFunctionBegin; 1037 PetscCall(TSReset_GLLE(ts)); 1038 if (ts->dm) { 1039 PetscCall(DMCoarsenHookRemove(ts->dm, DMCoarsenHook_TSGLLE, DMRestrictHook_TSGLLE, ts)); 1040 PetscCall(DMSubDomainHookRemove(ts->dm, DMSubDomainHook_TSGLLE, DMSubDomainRestrictHook_TSGLLE, ts)); 1041 } 1042 if (gl->adapt) PetscCall(TSGLLEAdaptDestroy(&gl->adapt)); 1043 if (gl->Destroy) PetscCall((*gl->Destroy)(gl)); 1044 PetscCall(PetscFree(ts->data)); 1045 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSGLLESetType_C", NULL)); 1046 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSGLLESetAcceptType_C", NULL)); 1047 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSGLLEGetAdapt_C", NULL)); 1048 PetscFunctionReturn(PETSC_SUCCESS); 1049 } 1050 1051 /* 1052 This defines the nonlinear equation that is to be solved with SNES 1053 g(x) = f(t,x,z+shift*x) = 0 1054 */ 1055 static PetscErrorCode SNESTSFormFunction_GLLE(SNES snes, Vec x, Vec f, TS ts) 1056 { 1057 TS_GLLE *gl = (TS_GLLE *)ts->data; 1058 Vec Z, Ydot; 1059 DM dm, dmsave; 1060 1061 PetscFunctionBegin; 1062 PetscCall(SNESGetDM(snes, &dm)); 1063 PetscCall(TSGLLEGetVecs(ts, dm, &Z, &Ydot)); 1064 PetscCall(VecWAXPY(Ydot, gl->scoeff / ts->time_step, x, Z)); 1065 dmsave = ts->dm; 1066 ts->dm = dm; 1067 PetscCall(TSComputeIFunction(ts, gl->stage_time, x, Ydot, f, PETSC_FALSE)); 1068 ts->dm = dmsave; 1069 PetscCall(TSGLLERestoreVecs(ts, dm, &Z, &Ydot)); 1070 PetscFunctionReturn(PETSC_SUCCESS); 1071 } 1072 1073 static PetscErrorCode SNESTSFormJacobian_GLLE(SNES snes, Vec x, Mat A, Mat B, TS ts) 1074 { 1075 TS_GLLE *gl = (TS_GLLE *)ts->data; 1076 Vec Z, Ydot; 1077 DM dm, dmsave; 1078 1079 PetscFunctionBegin; 1080 PetscCall(SNESGetDM(snes, &dm)); 1081 PetscCall(TSGLLEGetVecs(ts, dm, &Z, &Ydot)); 1082 dmsave = ts->dm; 1083 ts->dm = dm; 1084 /* gl->Xdot will have already been computed in SNESTSFormFunction_GLLE */ 1085 PetscCall(TSComputeIJacobian(ts, gl->stage_time, x, gl->Ydot[gl->stage], gl->scoeff / ts->time_step, A, B, PETSC_FALSE)); 1086 ts->dm = dmsave; 1087 PetscCall(TSGLLERestoreVecs(ts, dm, &Z, &Ydot)); 1088 PetscFunctionReturn(PETSC_SUCCESS); 1089 } 1090 1091 static PetscErrorCode TSSetUp_GLLE(TS ts) 1092 { 1093 TS_GLLE *gl = (TS_GLLE *)ts->data; 1094 PetscInt max_r, max_s; 1095 DM dm; 1096 1097 PetscFunctionBegin; 1098 gl->setupcalled = PETSC_TRUE; 1099 PetscCall(TSGLLEGetMaxSizes(ts, &max_r, &max_s)); 1100 PetscCall(VecDuplicateVecs(ts->vec_sol, max_r, &gl->X)); 1101 PetscCall(VecDuplicateVecs(ts->vec_sol, max_r, &gl->Xold)); 1102 PetscCall(VecDuplicateVecs(ts->vec_sol, max_s, &gl->Ydot)); 1103 PetscCall(VecDuplicateVecs(ts->vec_sol, 3, &gl->himom)); 1104 PetscCall(VecDuplicate(ts->vec_sol, &gl->W)); 1105 PetscCall(VecDuplicate(ts->vec_sol, &gl->Y)); 1106 PetscCall(VecDuplicate(ts->vec_sol, &gl->Z)); 1107 1108 /* Default acceptance tests and adaptivity */ 1109 if (!gl->Accept) PetscCall(TSGLLESetAcceptType(ts, TSGLLEACCEPT_ALWAYS)); 1110 if (!gl->adapt) PetscCall(TSGLLEGetAdapt(ts, &gl->adapt)); 1111 1112 if (gl->current_scheme < 0) { 1113 PetscInt i; 1114 for (i = 0;; i++) { 1115 if (gl->schemes[i]->p == gl->start_order) break; 1116 PetscCheck(i + 1 != gl->nschemes, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No schemes available with requested start order %" PetscInt_FMT, i); 1117 } 1118 gl->current_scheme = i; 1119 } 1120 PetscCall(TSGetDM(ts, &dm)); 1121 PetscCall(DMCoarsenHookAdd(dm, DMCoarsenHook_TSGLLE, DMRestrictHook_TSGLLE, ts)); 1122 PetscCall(DMSubDomainHookAdd(dm, DMSubDomainHook_TSGLLE, DMSubDomainRestrictHook_TSGLLE, ts)); 1123 PetscFunctionReturn(PETSC_SUCCESS); 1124 } 1125 /*------------------------------------------------------------*/ 1126 1127 static PetscErrorCode TSSetFromOptions_GLLE(TS ts, PetscOptionItems *PetscOptionsObject) 1128 { 1129 TS_GLLE *gl = (TS_GLLE *)ts->data; 1130 char tname[256] = TSGLLE_IRKS, completef[256] = "rescale-and-modify"; 1131 1132 PetscFunctionBegin; 1133 PetscOptionsHeadBegin(PetscOptionsObject, "General Linear ODE solver options"); 1134 { 1135 PetscBool flg; 1136 PetscCall(PetscOptionsFList("-ts_gl_type", "Type of GL method", "TSGLLESetType", TSGLLEList, gl->type_name[0] ? gl->type_name : tname, tname, sizeof(tname), &flg)); 1137 if (flg || !gl->type_name[0]) PetscCall(TSGLLESetType(ts, tname)); 1138 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)); 1139 PetscCall(PetscOptionsInt("-ts_gl_max_order", "Maximum order to try", "TSGLLESetMaxOrder", gl->max_order, &gl->max_order, NULL)); 1140 PetscCall(PetscOptionsInt("-ts_gl_min_order", "Minimum order to try", "TSGLLESetMinOrder", gl->min_order, &gl->min_order, NULL)); 1141 PetscCall(PetscOptionsInt("-ts_gl_start_order", "Initial order to try", "TSGLLESetMinOrder", gl->start_order, &gl->start_order, NULL)); 1142 PetscCall(PetscOptionsEnum("-ts_gl_error_direction", "Which direction to look when estimating error", "TSGLLESetErrorDirection", TSGLLEErrorDirections, (PetscEnum)gl->error_direction, (PetscEnum *)&gl->error_direction, NULL)); 1143 PetscCall(PetscOptionsBool("-ts_gl_extrapolate", "Extrapolate stage solution from previous solution (sometimes unstable)", "TSGLLESetExtrapolate", gl->extrapolate, &gl->extrapolate, NULL)); 1144 PetscCall(PetscOptionsReal("-ts_gl_atol", "Absolute tolerance", "TSGLLESetTolerances", gl->wrms_atol, &gl->wrms_atol, NULL)); 1145 PetscCall(PetscOptionsReal("-ts_gl_rtol", "Relative tolerance", "TSGLLESetTolerances", gl->wrms_rtol, &gl->wrms_rtol, NULL)); 1146 PetscCall(PetscOptionsString("-ts_gl_complete", "Method to use for completing the step", "none", completef, completef, sizeof(completef), &flg)); 1147 if (flg) { 1148 PetscBool match1, match2; 1149 PetscCall(PetscStrcmp(completef, "rescale", &match1)); 1150 PetscCall(PetscStrcmp(completef, "rescale-and-modify", &match2)); 1151 if (match1) gl->CompleteStep = TSGLLECompleteStep_Rescale; 1152 else if (match2) gl->CompleteStep = TSGLLECompleteStep_RescaleAndModify; 1153 else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "%s", completef); 1154 } 1155 { 1156 char type[256] = TSGLLEACCEPT_ALWAYS; 1157 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)); 1158 if (flg || !gl->accept_name[0]) PetscCall(TSGLLESetAcceptType(ts, type)); 1159 } 1160 { 1161 TSGLLEAdapt adapt; 1162 PetscCall(TSGLLEGetAdapt(ts, &adapt)); 1163 PetscCall(TSGLLEAdaptSetFromOptions(adapt, PetscOptionsObject)); 1164 } 1165 } 1166 PetscOptionsHeadEnd(); 1167 PetscFunctionReturn(PETSC_SUCCESS); 1168 } 1169 1170 static PetscErrorCode TSView_GLLE(TS ts, PetscViewer viewer) 1171 { 1172 TS_GLLE *gl = (TS_GLLE *)ts->data; 1173 PetscInt i; 1174 PetscBool iascii, details; 1175 1176 PetscFunctionBegin; 1177 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii)); 1178 if (iascii) { 1179 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)); 1180 PetscCall(PetscViewerASCIIPrintf(viewer, " Error estimation: %s\n", TSGLLEErrorDirections[gl->error_direction])); 1181 PetscCall(PetscViewerASCIIPrintf(viewer, " Extrapolation: %s\n", gl->extrapolate ? "yes" : "no")); 1182 PetscCall(PetscViewerASCIIPrintf(viewer, " Acceptance test: %s\n", gl->accept_name[0] ? gl->accept_name : "(not yet set)")); 1183 PetscCall(PetscViewerASCIIPushTab(viewer)); 1184 PetscCall(TSGLLEAdaptView(gl->adapt, viewer)); 1185 PetscCall(PetscViewerASCIIPopTab(viewer)); 1186 PetscCall(PetscViewerASCIIPrintf(viewer, " type: %s\n", gl->type_name[0] ? gl->type_name : "(not yet set)")); 1187 PetscCall(PetscViewerASCIIPrintf(viewer, "Schemes within family (%" PetscInt_FMT "):\n", gl->nschemes)); 1188 details = PETSC_FALSE; 1189 PetscCall(PetscOptionsGetBool(((PetscObject)ts)->options, ((PetscObject)ts)->prefix, "-ts_gl_view_detailed", &details, NULL)); 1190 PetscCall(PetscViewerASCIIPushTab(viewer)); 1191 for (i = 0; i < gl->nschemes; i++) PetscCall(TSGLLESchemeView(gl->schemes[i], details, viewer)); 1192 if (gl->View) PetscCall((*gl->View)(gl, viewer)); 1193 PetscCall(PetscViewerASCIIPopTab(viewer)); 1194 } 1195 PetscFunctionReturn(PETSC_SUCCESS); 1196 } 1197 1198 /*@C 1199 TSGLLERegister - adds a `TSGLLE` implementation 1200 1201 Not Collective 1202 1203 Input Parameters: 1204 + sname - name of user-defined general linear scheme 1205 - function - routine to create method context 1206 1207 Level: advanced 1208 1209 Note: 1210 `TSGLLERegister()` may be called multiple times to add several user-defined families. 1211 1212 Sample usage: 1213 .vb 1214 TSGLLERegister("my_scheme",MySchemeCreate); 1215 .ve 1216 1217 Then, your scheme can be chosen with the procedural interface via 1218 $ TSGLLESetType(ts,"my_scheme") 1219 or at runtime via the option 1220 $ -ts_gl_type my_scheme 1221 1222 .seealso: [](chapter_ts), `TSGLLE`, `TSGLLEType`, `TSGLLERegisterAll()` 1223 @*/ 1224 PetscErrorCode TSGLLERegister(const char sname[], PetscErrorCode (*function)(TS)) 1225 { 1226 PetscFunctionBegin; 1227 PetscCall(TSGLLEInitializePackage()); 1228 PetscCall(PetscFunctionListAdd(&TSGLLEList, sname, function)); 1229 PetscFunctionReturn(PETSC_SUCCESS); 1230 } 1231 1232 /*@C 1233 TSGLLEAcceptRegister - adds a `TSGLLE` acceptance scheme 1234 1235 Not Collective 1236 1237 Input Parameters: 1238 + sname - name of user-defined acceptance scheme 1239 - function - routine to create method context 1240 1241 Level: advanced 1242 1243 Note: 1244 `TSGLLEAcceptRegister()` may be called multiple times to add several user-defined families. 1245 1246 Sample usage: 1247 .vb 1248 TSGLLEAcceptRegister("my_scheme",MySchemeCreate); 1249 .ve 1250 1251 Then, your scheme can be chosen with the procedural interface via 1252 $ TSGLLESetAcceptType(ts,"my_scheme") 1253 or at runtime via the option 1254 $ -ts_gl_accept_type my_scheme 1255 1256 .seealso: [](chapter_ts), `TSGLLE`, `TSGLLEType`, `TSGLLERegisterAll()`, `TSGLLEAcceptFunction` 1257 @*/ 1258 PetscErrorCode TSGLLEAcceptRegister(const char sname[], TSGLLEAcceptFunction 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: [](chapter_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: [](chapter_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: [](chapter_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 1324 1325 These methods contain Runge-Kutta and multistep schemes as special cases. These special cases have some fundamental 1326 limitations. For example, diagonally implicit Runge-Kutta cannot have stage order greater than 1 which limits their 1327 applicability to very stiff systems. Meanwhile, multistep methods cannot be A-stable for order greater than 2 and BDF 1328 are not 0-stable for order greater than 6. GL methods can be A- and L-stable with arbitrarily high stage order and 1329 reliable error estimates for both 1 and 2 orders higher to facilitate adaptive step sizes and adaptive order schemes. 1330 All this is possible while preserving a singly diagonally implicit structure. 1331 1332 Options Database Keys: 1333 + -ts_gl_type <type> - the class of general linear method (irks) 1334 . -ts_gl_rtol <tol> - relative error 1335 . -ts_gl_atol <tol> - absolute error 1336 . -ts_gl_min_order <p> - minimum order method to consider (default=1) 1337 . -ts_gl_max_order <p> - maximum order method to consider (default=3) 1338 . -ts_gl_start_order <p> - order of starting method (default=1) 1339 . -ts_gl_complete <method> - method to use for completing the step (rescale-and-modify or rescale) 1340 - -ts_adapt_type <method> - adaptive controller to use (none step both) 1341 1342 Level: beginner 1343 1344 Notes: 1345 This integrator can be applied to DAE. 1346 1347 Diagonally implicit general linear (DIGL) methods are a generalization of diagonally implicit Runge-Kutta (DIRK). 1348 They are represented by the tableau 1349 1350 .vb 1351 A | U 1352 ------- 1353 B | V 1354 .ve 1355 1356 combined with a vector c of abscissa. "Diagonally implicit" means that A is lower triangular. 1357 A step of the general method reads 1358 1359 .vb 1360 [ Y ] = [A U] [ Y' ] 1361 [X^k] = [B V] [X^{k-1}] 1362 .ve 1363 1364 where Y is the multivector of stage values, Y' is the multivector of stage derivatives, X^k is the Nordsieck vector of 1365 the solution at step k. The Nordsieck vector consists of the first r moments of the solution, given by 1366 1367 .vb 1368 X = [x_0,x_1,...,x_{r-1}] = [x, h x', h^2 x'', ..., h^{r-1} x^{(r-1)} ] 1369 .ve 1370 1371 If A is lower triangular, we can solve the stages (Y,Y') sequentially 1372 1373 .vb 1374 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} 1375 .ve 1376 1377 and then construct the pieces to carry to the next step 1378 1379 .vb 1380 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} 1381 .ve 1382 1383 Note that when the equations are cast in implicit form, we are using the stage equation to define y'_i 1384 in terms of y_i and known stuff (y_j for j<i and x_j for all j). 1385 1386 Error estimation 1387 1388 At present, the most attractive GL methods for stiff problems are singly diagonally implicit schemes which posses 1389 Inherent Runge-Kutta Stability (`TSIRKS`). These methods have r=s, the number of items passed between steps is equal to 1390 the number of stages. The order and stage-order are one less than the number of stages. We use the error estimates 1391 in the 2007 paper which provide the following estimates 1392 1393 .vb 1394 h^{p+1} X^{(p+1)} = phi_0^T Y' + [0 psi_0^T] Xold 1395 h^{p+2} X^{(p+2)} = phi_1^T Y' + [0 psi_1^T] Xold 1396 h^{p+2} (dx'/dx) X^{(p+1)} = phi_2^T Y' + [0 psi_2^T] Xold 1397 .ve 1398 1399 These estimates are accurate to O(h^{p+3}). 1400 1401 Changing the step size 1402 1403 We use the generalized "rescale and modify" scheme, see equation (4.5) of the 2007 paper. 1404 1405 References: 1406 + * - John Butcher and Z. Jackieweicz and W. Wright, On error propagation in general linear methods for 1407 ordinary differential equations, Journal of Complexity, Vol 23, 2007. 1408 - * - John Butcher, Numerical methods for ordinary differential equations, second edition, Wiley, 2009. 1409 1410 .seealso: [](chapter_ts), `TSCreate()`, `TS`, `TSSetType()`, `TSType` 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(PETSC_SUCCESS); 1447 } 1448