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 CHKERRQ(DMGetNamedGlobalVector(dm,"TSGLLE_Z",Z)); 41 } else *Z = gl->Z; 42 } 43 if (Ydotstage) { 44 if (dm && dm != ts->dm) { 45 CHKERRQ(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 CHKERRQ(DMRestoreNamedGlobalVector(dm,"TSGLLE_Z",Z)); 57 } 58 } 59 if (Ydotstage) { 60 61 if (dm && dm != ts->dm) { 62 CHKERRQ(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 CHKERRQ(TSGLLEGetVecs(ts,fine,NULL,&Ydot)); 81 CHKERRQ(TSGLLEGetVecs(ts,coarse,NULL,&Ydot_c)); 82 CHKERRQ(MatRestrict(restrct,Ydot,Ydot_c)); 83 CHKERRQ(VecPointwiseMult(Ydot_c,rscale,Ydot_c)); 84 CHKERRQ(TSGLLERestoreVecs(ts,fine,NULL,&Ydot)); 85 CHKERRQ(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 CHKERRQ(TSGLLEGetVecs(ts,dm,NULL,&Ydot)); 102 CHKERRQ(TSGLLEGetVecs(ts,subdm,NULL,&Ydot_s)); 103 104 CHKERRQ(VecScatterBegin(gscat,Ydot,Ydot_s,INSERT_VALUES,SCATTER_FORWARD)); 105 CHKERRQ(VecScatterEnd(gscat,Ydot,Ydot_s,INSERT_VALUES,SCATTER_FORWARD)); 106 107 CHKERRQ(TSGLLERestoreVecs(ts,dm,NULL,&Ydot)); 108 CHKERRQ(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 PetscCheckFalse(p < 1,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Scheme order must be positive"); 120 PetscCheckFalse(r < 1,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"At least one item must be carried between steps"); 121 PetscCheckFalse(s < 1,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"At least one stage is required"); 122 PetscValidPointer(inscheme,10); 123 *inscheme = NULL; 124 CHKERRQ(PetscNew(&scheme)); 125 scheme->p = p; 126 scheme->q = q; 127 scheme->r = r; 128 scheme->s = s; 129 130 CHKERRQ(PetscMalloc5(s,&scheme->c,s*s,&scheme->a,r*s,&scheme->b,r*s,&scheme->u,r*r,&scheme->v)); 131 CHKERRQ(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 CHKERRQ(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 CHKERRQ(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 CHKERRQ(PetscBLASIntCast(r-1,&m)); 156 CHKERRQ(PetscBLASIntCast(r,&n)); 157 PetscStackCallBLAS("LAPACKgesv",LAPACKgesv_(&m,&one,ImV,&n,ipiv,scheme->alpha+1,&n,&info)); 158 PetscCheckFalse(info < 0,PETSC_COMM_SELF,PETSC_ERR_LIB,"Bad argument to GESV"); 159 PetscCheckFalse(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 PetscCheckFalse(info < 0,PETSC_COMM_SELF,PETSC_ERR_LIB,"Bad argument to GETRS"); 168 PetscCheckFalse(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 PetscCheckFalse(info < 0,PETSC_COMM_SELF,PETSC_ERR_LIB,"Bad argument to GETRS"); 193 PetscCheckFalse(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 * % Determine 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 CHKERRQ(PetscBLASIntCast(s,&n)); 237 CHKERRQ(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 PetscCheckFalse(info < 0,PETSC_COMM_SELF,PETSC_ERR_LIB,"Bad argument to GELSS"); 247 PetscCheckFalse(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 CHKERRQ(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 CHKERRQ(PetscFree5(sc->c,sc->a,sc->b,sc->u,sc->v)); 287 CHKERRQ(PetscFree6(sc->alpha,sc->beta,sc->gamma,sc->phi,sc->psi,sc->stage_error)); 288 CHKERRQ(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]) CHKERRQ(TSGLLESchemeDestroy(gl->schemes[i])); 299 } 300 CHKERRQ(PetscFree(gl->schemes)); 301 gl->nschemes = 0; 302 CHKERRQ(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 CHKERRQ(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii)); 313 if (iascii) { 314 CHKERRQ(PetscViewerASCIIPrintf(viewer,"%30s = [",name)); 315 for (i=0; i<m; i++) { 316 if (i) CHKERRQ(PetscViewerASCIIPrintf(viewer,"%30s [","")); 317 CHKERRQ(PetscViewerASCIIUseTabs(viewer,PETSC_FALSE)); 318 for (j=0; j<n; j++) { 319 CHKERRQ(PetscViewerASCIIPrintf(viewer," %12.8g",PetscRealPart(a[i*n+j]))); 320 } 321 CHKERRQ(PetscViewerASCIIPrintf(viewer,"]\n")); 322 CHKERRQ(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 CHKERRQ(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii)); 334 if (iascii) { 335 CHKERRQ(PetscViewerASCIIPrintf(viewer,"GL scheme p,q,r,s = %d,%d,%d,%d\n",sc->p,sc->q,sc->r,sc->s)); 336 CHKERRQ(PetscViewerASCIIPushTab(viewer)); 337 CHKERRQ(PetscViewerASCIIPrintf(viewer,"Stiffly accurate: %s, FSAL: %s\n",sc->stiffly_accurate ? "yes" : "no",sc->fsal ? "yes" : "no")); 338 CHKERRQ(PetscViewerASCIIPrintf(viewer,"Leading error constants: %10.3e %10.3e %10.3e\n", 339 PetscRealPart(sc->alpha[0]),PetscRealPart(sc->beta[0]),PetscRealPart(sc->gamma[0]))); 340 CHKERRQ(TSGLLEViewTable_Private(viewer,1,sc->s,sc->c,"Abscissas c")); 341 if (view_details) { 342 CHKERRQ(TSGLLEViewTable_Private(viewer,sc->s,sc->s,sc->a,"A")); 343 CHKERRQ(TSGLLEViewTable_Private(viewer,sc->r,sc->s,sc->b,"B")); 344 CHKERRQ(TSGLLEViewTable_Private(viewer,sc->s,sc->r,sc->u,"U")); 345 CHKERRQ(TSGLLEViewTable_Private(viewer,sc->r,sc->r,sc->v,"V")); 346 347 CHKERRQ(TSGLLEViewTable_Private(viewer,3,sc->s,sc->phi,"Error estimate phi")); 348 CHKERRQ(TSGLLEViewTable_Private(viewer,3,sc->r,sc->psi,"Error estimate psi")); 349 CHKERRQ(TSGLLEViewTable_Private(viewer,1,sc->r,sc->alpha,"Modify alpha")); 350 CHKERRQ(TSGLLEViewTable_Private(viewer,1,sc->r,sc->beta,"Modify beta")); 351 CHKERRQ(TSGLLEViewTable_Private(viewer,1,sc->r,sc->gamma,"Modify gamma")); 352 CHKERRQ(TSGLLEViewTable_Private(viewer,1,sc->s,sc->stage_error,"Stage error xi")); 353 } 354 CHKERRQ(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 PetscCheckFalse(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 CHKERRQ(VecZeroEntries(hm[i])); 371 CHKERRQ(VecMAXPY(hm[i],sc->s,phih,Ydot)); 372 CHKERRQ(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 CHKERRQ(VecZeroEntries(X[i])); 388 for (j=0; j<s; j++) brow[j] = h*sc->b[i*s+j]; 389 CHKERRQ(VecMAXPY(X[i],s,brow,Ydot)); 390 for (j=0; j<r; j++) vrow[j] = sc->v[i*r+j]; 391 CHKERRQ(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 CHKERRQ(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 CHKERRQ(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 CHKERRQ(VecMAXPY(X[i],r,vrow,Xold)); 424 } 425 if (r < next_sc->r) { 426 PetscCheckFalse(r+1 != next_sc->r,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Cannot accommodate jump in r greater than 1"); 427 CHKERRQ(VecZeroEntries(X[r])); 428 for (j=0; j<s; j++) brow[j] = h*PetscPowRealInt(ratio,p+1)*sc->phi[0*s+j]; 429 CHKERRQ(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 CHKERRQ(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 CHKERRQ(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 CHKERRQ(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 CHKERRQ(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 CHKERRQ(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 CHKERRQ(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 CHKERRQ(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 CHKERRQ(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 CHKERRQ(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 CHKERRQ(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 CHKERRQ(VecGetArray(gl->X[0],&x)); 680 CHKERRQ(VecGetArray(gl->W,&w)); 681 CHKERRQ(VecGetLocalSize(gl->W,&n)); 682 for (i=0; i<n; i++) w[i] = 1./(gl->wrms_atol + gl->wrms_rtol*PetscAbsScalar(x[i])); 683 CHKERRQ(VecRestoreArray(gl->X[0],&x)); 684 CHKERRQ(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 CHKERRQ(VecGetArray(X,&x)); 697 CHKERRQ(VecGetArray(gl->W,&w)); 698 CHKERRQ(VecGetLocalSize(gl->W,&n)); 699 for (i=0; i<n; i++) sum += PetscAbsScalar(PetscSqr(x[i]*w[i])); 700 CHKERRQ(VecRestoreArray(X,&x)); 701 CHKERRQ(VecRestoreArray(gl->W,&w)); 702 CHKERRMPI(MPIU_Allreduce(&sum,&gsum,1,MPIU_REAL,MPIU_SUM,PetscObjectComm((PetscObject)ts))); 703 CHKERRQ(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 CHKERRQ(PetscStrcmp(gl->type_name,type,&same)); 717 if (same) PetscFunctionReturn(0); 718 CHKERRQ((*gl->Destroy)(gl)); 719 } 720 721 CHKERRQ(PetscFunctionListFind(TSGLLEList,type,&r)); 722 PetscCheck(r,PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown TSGLLE type \"%s\" given",type); 723 CHKERRQ((*r)(ts)); 724 CHKERRQ(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 CHKERRQ(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 CHKERRQ(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 CHKERRQ(TSGLLEAdaptCreate(PetscObjectComm((PetscObject)ts),&gl->adapt)); 748 CHKERRQ(PetscObjectIncrementTabLevel((PetscObject)gl->adapt,(PetscObject)ts,1)); 749 CHKERRQ(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 PetscCheckFalse(cur < 0 || gl->nschemes <= cur,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Current scheme not found in scheme list"); 779 CHKERRQ(TSGLLEAdaptChoose(gl->adapt,n,orders,errors,costs,cur,h,tleft,&next_sc,next_h,finish)); 780 *next_scheme = candidates[next_sc]; 781 CHKERRQ(PetscInfo(ts,"Adapt chose scheme %d (%d,%d,%d,%d) with step size %6.2e, finish=%d\n",*next_scheme,gl->schemes[*next_scheme]->p,gl->schemes[*next_scheme]->q,gl->schemes[*next_scheme]->r,gl->schemes[*next_scheme]->s,*next_h,*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 CHKERRQ(TSMonitor(ts,ts->steps,ts->ptime,ts->vec_sol)); 804 805 CHKERRQ(TSGLLEGetMaxSizes(ts,&max_r,&max_s)); 806 CHKERRQ(VecCopy(ts->vec_sol,gl->X[0])); 807 for (i=1; i<max_r; i++) { 808 CHKERRQ(VecZeroEntries(gl->X[i])); 809 } 810 CHKERRQ(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 CHKERRQ(VecCopy(ts->vec_sol,gl->Z)); 819 CHKERRQ(VecCopy(ts->vec_sol,gl->Y)); 820 CHKERRQ(SNESSolve(ts->snes,NULL,gl->Y)); 821 CHKERRQ(SNESGetIterationNumber(ts->snes,&its)); 822 CHKERRQ(SNESGetLinearSolveIterations(ts->snes,&lits)); 823 CHKERRQ(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 CHKERRQ(PetscInfo(ts,"Step=%D, nonlinear solve solve failures %D greater than current TS allowed, stopping solve\n",ts->steps,ts->num_snes_failures)); 829 PetscFunctionReturn(0); 830 } 831 } 832 833 PetscCheckFalse(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 CHKERRQ(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 CHKERRQ(VecZeroEntries(gl->Z)); 873 for (j=0; j<r; j++) { 874 CHKERRQ(VecAXPY(gl->Z,-shift*u[i*r+j],X[j])); 875 } 876 for (j=0; j<i; j++) { 877 CHKERRQ(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 CHKERRQ(VecWAXPY(Y,c[i]*h,X[1],X[0])); 886 } else { 887 /* Linear extrapolation from the last stage */ 888 CHKERRQ(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 CHKERRQ(VecCopy(X[0],Y)); 892 } 893 894 /* Solve this stage (Ydot[i] is computed during function evaluation) */ 895 CHKERRQ(SNESSolve(ts->snes,NULL,Y)); 896 CHKERRQ(SNESGetIterationNumber(ts->snes,&its)); 897 CHKERRQ(SNESGetLinearSolveIterations(ts->snes,&lits)); 898 CHKERRQ(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 CHKERRQ(PetscInfo(ts,"Step=%D, nonlinear solve solve failures %D 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 CHKERRQ((*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 CHKERRQ(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 CHKERRQ((*gl->Accept)(ts,ts->max_time-gl->stage_time,h,enorm,&accept)); 918 if (accept) goto accepted; 919 rejections++; 920 CHKERRQ(PetscInfo(ts,"Step %D (t=%g) not accepted, rejections=%D\n",k,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 CHKERRQ(VecScale(X[i],PetscPowRealInt(0.5,i))); 933 } 934 } 935 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_CONV_FAILED,"Time step %D (t=%g) not accepted after %D failures",k,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 CHKERRQ(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 CHKERRQ(PetscInfo(ts,"Last moment norm %10.2e, estimated error norms %10.2e %10.2e %10.2e\n",hmnorm[0],enorm[0],enorm[1],enorm[2])); 943 if (!final_step) { 944 CHKERRQ(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 CHKERRQ((*gl->CompleteStep)(scheme,h,gl->schemes[next_scheme],next_h,Ydot,gl->Xold,gl->X)); 956 957 CHKERRQ(TSGLLEUpdateWRMS(ts)); 958 959 /* Post the solution for the user, we could avoid this copy with a small bit of cleverness */ 960 CHKERRQ(VecCopy(gl->X[0],ts->vec_sol)); 961 ts->ptime += h; 962 ts->steps++; 963 964 CHKERRQ(TSPostEvaluate(ts)); 965 CHKERRQ(TSPostStep(ts)); 966 CHKERRQ(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 CHKERRQ(TSGLLEGetMaxSizes(ts,&max_r,&max_s)); 984 CHKERRQ(VecDestroyVecs(max_r,&gl->Xold)); 985 CHKERRQ(VecDestroyVecs(max_r,&gl->X)); 986 CHKERRQ(VecDestroyVecs(max_s,&gl->Ydot)); 987 CHKERRQ(VecDestroyVecs(3,&gl->himom)); 988 CHKERRQ(VecDestroy(&gl->W)); 989 CHKERRQ(VecDestroy(&gl->Y)); 990 CHKERRQ(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 CHKERRQ(TSReset_GLLE(ts)); 1002 if (ts->dm) { 1003 CHKERRQ(DMCoarsenHookRemove(ts->dm,DMCoarsenHook_TSGLLE,DMRestrictHook_TSGLLE,ts)); 1004 CHKERRQ(DMSubDomainHookRemove(ts->dm,DMSubDomainHook_TSGLLE,DMSubDomainRestrictHook_TSGLLE,ts)); 1005 } 1006 if (gl->adapt) CHKERRQ(TSGLLEAdaptDestroy(&gl->adapt)); 1007 if (gl->Destroy) CHKERRQ((*gl->Destroy)(gl)); 1008 CHKERRQ(PetscFree(ts->data)); 1009 CHKERRQ(PetscObjectComposeFunction((PetscObject)ts,"TSGLLESetType_C",NULL)); 1010 CHKERRQ(PetscObjectComposeFunction((PetscObject)ts,"TSGLLESetAcceptType_C",NULL)); 1011 CHKERRQ(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 CHKERRQ(SNESGetDM(snes,&dm)); 1027 CHKERRQ(TSGLLEGetVecs(ts,dm,&Z,&Ydot)); 1028 CHKERRQ(VecWAXPY(Ydot,gl->scoeff/ts->time_step,x,Z)); 1029 dmsave = ts->dm; 1030 ts->dm = dm; 1031 CHKERRQ(TSComputeIFunction(ts,gl->stage_time,x,Ydot,f,PETSC_FALSE)); 1032 ts->dm = dmsave; 1033 CHKERRQ(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 CHKERRQ(SNESGetDM(snes,&dm)); 1045 CHKERRQ(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 CHKERRQ(TSComputeIJacobian(ts,gl->stage_time,x,gl->Ydot[gl->stage],gl->scoeff/ts->time_step,A,B,PETSC_FALSE)); 1050 ts->dm = dmsave; 1051 CHKERRQ(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 CHKERRQ(TSGLLEGetMaxSizes(ts,&max_r,&max_s)); 1064 CHKERRQ(VecDuplicateVecs(ts->vec_sol,max_r,&gl->X)); 1065 CHKERRQ(VecDuplicateVecs(ts->vec_sol,max_r,&gl->Xold)); 1066 CHKERRQ(VecDuplicateVecs(ts->vec_sol,max_s,&gl->Ydot)); 1067 CHKERRQ(VecDuplicateVecs(ts->vec_sol,3,&gl->himom)); 1068 CHKERRQ(VecDuplicate(ts->vec_sol,&gl->W)); 1069 CHKERRQ(VecDuplicate(ts->vec_sol,&gl->Y)); 1070 CHKERRQ(VecDuplicate(ts->vec_sol,&gl->Z)); 1071 1072 /* Default acceptance tests and adaptivity */ 1073 if (!gl->Accept) CHKERRQ(TSGLLESetAcceptType(ts,TSGLLEACCEPT_ALWAYS)); 1074 if (!gl->adapt) CHKERRQ(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 PetscCheckFalse(i+1 == gl->nschemes,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"No schemes available with requested start order %d",i); 1081 } 1082 gl->current_scheme = i; 1083 } 1084 CHKERRQ(TSGetDM(ts,&dm)); 1085 CHKERRQ(DMCoarsenHookAdd(dm,DMCoarsenHook_TSGLLE,DMRestrictHook_TSGLLE,ts)); 1086 CHKERRQ(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 CHKERRQ(PetscOptionsHead(PetscOptionsObject,"General Linear ODE solver options")); 1098 { 1099 PetscBool flg; 1100 CHKERRQ(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 CHKERRQ(TSGLLESetType(ts,tname)); 1103 } 1104 CHKERRQ(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 CHKERRQ(PetscOptionsInt("-ts_gl_max_order","Maximum order to try","TSGLLESetMaxOrder",gl->max_order,&gl->max_order,NULL)); 1106 CHKERRQ(PetscOptionsInt("-ts_gl_min_order","Minimum order to try","TSGLLESetMinOrder",gl->min_order,&gl->min_order,NULL)); 1107 CHKERRQ(PetscOptionsInt("-ts_gl_start_order","Initial order to try","TSGLLESetMinOrder",gl->start_order,&gl->start_order,NULL)); 1108 CHKERRQ(PetscOptionsEnum("-ts_gl_error_direction","Which direction to look when estimating error","TSGLLESetErrorDirection",TSGLLEErrorDirections,(PetscEnum)gl->error_direction,(PetscEnum*)&gl->error_direction,NULL)); 1109 CHKERRQ(PetscOptionsBool("-ts_gl_extrapolate","Extrapolate stage solution from previous solution (sometimes unstable)","TSGLLESetExtrapolate",gl->extrapolate,&gl->extrapolate,NULL)); 1110 CHKERRQ(PetscOptionsReal("-ts_gl_atol","Absolute tolerance","TSGLLESetTolerances",gl->wrms_atol,&gl->wrms_atol,NULL)); 1111 CHKERRQ(PetscOptionsReal("-ts_gl_rtol","Relative tolerance","TSGLLESetTolerances",gl->wrms_rtol,&gl->wrms_rtol,NULL)); 1112 CHKERRQ(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 CHKERRQ(PetscStrcmp(completef,"rescale",&match1)); 1116 CHKERRQ(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 CHKERRQ(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 CHKERRQ(TSGLLESetAcceptType(ts,type)); 1126 } 1127 } 1128 { 1129 TSGLLEAdapt adapt; 1130 CHKERRQ(TSGLLEGetAdapt(ts,&adapt)); 1131 CHKERRQ(TSGLLEAdaptSetFromOptions(PetscOptionsObject,adapt)); 1132 } 1133 } 1134 CHKERRQ(PetscOptionsTail()); 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 CHKERRQ(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii)); 1146 if (iascii) { 1147 CHKERRQ(PetscViewerASCIIPrintf(viewer," min order %D, max order %D, current order %D\n",gl->min_order,gl->max_order,gl->schemes[gl->current_scheme]->p)); 1148 CHKERRQ(PetscViewerASCIIPrintf(viewer," Error estimation: %s\n",TSGLLEErrorDirections[gl->error_direction])); 1149 CHKERRQ(PetscViewerASCIIPrintf(viewer," Extrapolation: %s\n",gl->extrapolate ? "yes" : "no")); 1150 CHKERRQ(PetscViewerASCIIPrintf(viewer," Acceptance test: %s\n",gl->accept_name[0] ? gl->accept_name : "(not yet set)")); 1151 CHKERRQ(PetscViewerASCIIPushTab(viewer)); 1152 CHKERRQ(TSGLLEAdaptView(gl->adapt,viewer)); 1153 CHKERRQ(PetscViewerASCIIPopTab(viewer)); 1154 CHKERRQ(PetscViewerASCIIPrintf(viewer," type: %s\n",gl->type_name[0] ? gl->type_name : "(not yet set)")); 1155 CHKERRQ(PetscViewerASCIIPrintf(viewer,"Schemes within family (%d):\n",gl->nschemes)); 1156 details = PETSC_FALSE; 1157 CHKERRQ(PetscOptionsGetBool(((PetscObject)ts)->options,((PetscObject)ts)->prefix,"-ts_gl_view_detailed",&details,NULL)); 1158 CHKERRQ(PetscViewerASCIIPushTab(viewer)); 1159 for (i=0; i<gl->nschemes; i++) { 1160 CHKERRQ(TSGLLESchemeView(gl->schemes[i],details,viewer)); 1161 } 1162 if (gl->View) { 1163 CHKERRQ((*gl->View)(gl,viewer)); 1164 } 1165 CHKERRQ(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 CHKERRQ(TSGLLEInitializePackage()); 1200 CHKERRQ(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 CHKERRQ(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 CHKERRQ(TSGLLERegister(TSGLLE_IRKS, TSGLLECreate_IRKS)); 1253 CHKERRQ(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 CHKERRQ(TSGLLERegisterAll()); 1271 CHKERRQ(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 CHKERRQ(PetscFunctionListDestroy(&TSGLLEList)); 1287 CHKERRQ(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 CHKERRQ(TSGLLEInitializePackage()); 1391 1392 CHKERRQ(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 CHKERRQ(PetscObjectComposeFunction((PetscObject)ts,"TSGLLESetType_C", &TSGLLESetType_GLLE)); 1417 CHKERRQ(PetscObjectComposeFunction((PetscObject)ts,"TSGLLESetAcceptType_C",&TSGLLESetAcceptType_GLLE)); 1418 CHKERRQ(PetscObjectComposeFunction((PetscObject)ts,"TSGLLEGetAdapt_C", &TSGLLEGetAdapt_GLLE)); 1419 PetscFunctionReturn(0); 1420 } 1421