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