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