#include <../src/ts/impls/implicit/glle/glle.h> /*I "petscts.h" I*/ #include #include static const char *TSGLLEErrorDirections[] = {"FORWARD","BACKWARD","TSGLLEErrorDirection","TSGLLEERROR_",0}; static PetscFunctionList TSGLLEList; static PetscFunctionList TSGLLEAcceptList; static PetscBool TSGLLEPackageInitialized; static PetscBool TSGLLERegisterAllCalled; /* This function is pure */ static PetscScalar Factorial(PetscInt n) { PetscInt i; if (n < 12) { /* Can compute with 32-bit integers */ PetscInt f = 1; for (i=2; i<=n; i++) f *= i; return (PetscScalar)f; } else { PetscScalar f = 1.; for (i=2; i<=n; i++) f *= (PetscScalar)i; return f; } } /* This function is pure */ static PetscScalar CPowF(PetscScalar c,PetscInt p) { return PetscPowRealInt(PetscRealPart(c),p)/Factorial(p); } #undef __FUNCT__ #define __FUNCT__ "TSGLLEGetVecs" static PetscErrorCode TSGLLEGetVecs(TS ts,DM dm,Vec *Z,Vec *Ydotstage) { TS_GLLE *gl = (TS_GLLE*)ts->data; PetscErrorCode ierr; PetscFunctionBegin; if (Z) { if (dm && dm != ts->dm) { ierr = DMGetNamedGlobalVector(dm,"TSGLLE_Z",Z);CHKERRQ(ierr); } else *Z = gl->Z; } if (Ydotstage) { if (dm && dm != ts->dm) { ierr = DMGetNamedGlobalVector(dm,"TSGLLE_Ydot",Ydotstage);CHKERRQ(ierr); } else *Ydotstage = gl->Ydot[gl->stage]; } PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "TSGLLERestoreVecs" static PetscErrorCode TSGLLERestoreVecs(TS ts,DM dm,Vec *Z,Vec *Ydotstage) { PetscErrorCode ierr; PetscFunctionBegin; if (Z) { if (dm && dm != ts->dm) { ierr = DMRestoreNamedGlobalVector(dm,"TSGLLE_Z",Z);CHKERRQ(ierr); } } if (Ydotstage) { if (dm && dm != ts->dm) { ierr = DMRestoreNamedGlobalVector(dm,"TSGLLE_Ydot",Ydotstage);CHKERRQ(ierr); } } PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "DMCoarsenHook_TSGLLE" static PetscErrorCode DMCoarsenHook_TSGLLE(DM fine,DM coarse,void *ctx) { PetscFunctionBegin; PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "DMRestrictHook_TSGLLE" static PetscErrorCode DMRestrictHook_TSGLLE(DM fine,Mat restrct,Vec rscale,Mat inject,DM coarse,void *ctx) { TS ts = (TS)ctx; PetscErrorCode ierr; Vec Ydot,Ydot_c; PetscFunctionBegin; ierr = TSGLLEGetVecs(ts,fine,NULL,&Ydot);CHKERRQ(ierr); ierr = TSGLLEGetVecs(ts,coarse,NULL,&Ydot_c);CHKERRQ(ierr); ierr = MatRestrict(restrct,Ydot,Ydot_c);CHKERRQ(ierr); ierr = VecPointwiseMult(Ydot_c,rscale,Ydot_c);CHKERRQ(ierr); ierr = TSGLLERestoreVecs(ts,fine,NULL,&Ydot);CHKERRQ(ierr); ierr = TSGLLERestoreVecs(ts,coarse,NULL,&Ydot_c);CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "DMSubDomainHook_TSGLLE" static PetscErrorCode DMSubDomainHook_TSGLLE(DM dm,DM subdm,void *ctx) { PetscFunctionBegin; PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "DMSubDomainRestrictHook_TSGLLE" static PetscErrorCode DMSubDomainRestrictHook_TSGLLE(DM dm,VecScatter gscat, VecScatter lscat,DM subdm,void *ctx) { TS ts = (TS)ctx; PetscErrorCode ierr; Vec Ydot,Ydot_s; PetscFunctionBegin; ierr = TSGLLEGetVecs(ts,dm,NULL,&Ydot);CHKERRQ(ierr); ierr = TSGLLEGetVecs(ts,subdm,NULL,&Ydot_s);CHKERRQ(ierr); ierr = VecScatterBegin(gscat,Ydot,Ydot_s,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); ierr = VecScatterEnd(gscat,Ydot,Ydot_s,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); ierr = TSGLLERestoreVecs(ts,dm,NULL,&Ydot);CHKERRQ(ierr); ierr = TSGLLERestoreVecs(ts,subdm,NULL,&Ydot_s);CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "TSGLLESchemeCreate" static PetscErrorCode TSGLLESchemeCreate(PetscInt p,PetscInt q,PetscInt r,PetscInt s,const PetscScalar *c, const PetscScalar *a,const PetscScalar *b,const PetscScalar *u,const PetscScalar *v,TSGLLEScheme *inscheme) { TSGLLEScheme scheme; PetscInt j; PetscErrorCode ierr; PetscFunctionBegin; if (p < 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Scheme order must be positive"); if (r < 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"At least one item must be carried between steps"); if (s < 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"At least one stage is required"); PetscValidPointer(inscheme,4); *inscheme = 0; ierr = PetscNew(&scheme);CHKERRQ(ierr); scheme->p = p; scheme->q = q; scheme->r = r; scheme->s = s; ierr = PetscMalloc5(s,&scheme->c,s*s,&scheme->a,r*s,&scheme->b,r*s,&scheme->u,r*r,&scheme->v);CHKERRQ(ierr); ierr = PetscMemcpy(scheme->c,c,s*sizeof(PetscScalar));CHKERRQ(ierr); for (j=0; ja[j] = (PetscAbsScalar(a[j]) < 1e-12) ? 0 : a[j]; for (j=0; jb[j] = (PetscAbsScalar(b[j]) < 1e-12) ? 0 : b[j]; for (j=0; ju[j] = (PetscAbsScalar(u[j]) < 1e-12) ? 0 : u[j]; for (j=0; jv[j] = (PetscAbsScalar(v[j]) < 1e-12) ? 0 : v[j]; ierr = PetscMalloc6(r,&scheme->alpha,r,&scheme->beta,r,&scheme->gamma,3*s,&scheme->phi,3*r,&scheme->psi,r,&scheme->stage_error);CHKERRQ(ierr); { PetscInt i,j,k,ss=s+2; PetscBLASInt m,n,one=1,*ipiv,lwork=4*((s+3)*3+3),info,ldb; PetscReal rcond,*sing,*workreal; PetscScalar *ImV,*H,*bmat,*workscalar,*c=scheme->c,*a=scheme->a,*b=scheme->b,*u=scheme->u,*v=scheme->v; #if !defined(PETSC_MISSING_LAPACK_GELSS) PetscBLASInt rank; #endif ierr = PetscMalloc7(PetscSqr(r),&ImV,3*s,&H,3*ss,&bmat,lwork,&workscalar,5*(3+r),&workreal,r+s,&sing,r+s,&ipiv);CHKERRQ(ierr); /* column-major input */ for (i=0; ialpha[i] = 1./Factorial(p+1-i); for (j=0; jalpha[i] -= b[i*s+j]*CPowF(c[j],p); } ierr = PetscBLASIntCast(r-1,&m);CHKERRQ(ierr); ierr = PetscBLASIntCast(r,&n);CHKERRQ(ierr); PetscStackCallBLAS("LAPACKgesv",LAPACKgesv_(&m,&one,ImV,&n,ipiv,scheme->alpha+1,&n,&info)); if (info < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Bad argument to GESV"); if (info > 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_MAT_LU_ZRPVT,"Bad LU factorization"); /* Build right hand side for beta (tp1 - glm.B(2:end,:)*(glm.c.^(p+1)./factorial(p+1)) - e.alpha) */ for (i=1; ibeta[i] = 1./Factorial(p+2-i) - scheme->alpha[i]; for (j=0; jbeta[i] -= b[i*s+j]*CPowF(c[j],p+1); } PetscStackCallBLAS("LAPACKgetrs",LAPACKgetrs_("No transpose",&m,&one,ImV,&n,ipiv,scheme->beta+1,&n,&info)); if (info < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Bad argument to GETRS"); if (info > 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Should not happen"); /* Build stage_error vector xi = glm.c.^(p+1)/factorial(p+1) - glm.A*glm.c.^p/factorial(p) + glm.U(:,2:end)*e.alpha; */ for (i=0; istage_error[i] = CPowF(c[i],p+1); for (j=0; jstage_error[i] -= a[i*s+j]*CPowF(c[j],p); for (j=1; jstage_error[i] += u[i*r+j]*scheme->alpha[j]; } /* alpha[0] (epsilon in B,J,W 2007) epsilon = 1/factorial(p+1) - B(1,:)*c.^p/factorial(p) + V(1,2:end)*e.alpha; */ scheme->alpha[0] = 1./Factorial(p+1); for (j=0; jalpha[0] -= b[0*s+j]*CPowF(c[j],p); for (j=1; jalpha[0] += v[0*r+j]*scheme->alpha[j]; /* right hand side for gamma (glm.B(2:end,:)*e.xi - e.epsilon*eye(s-1,1)) */ for (i=1; igamma[i] = (i==1 ? -1. : 0)*scheme->alpha[0]; for (j=0; jgamma[i] += b[i*s+j]*scheme->stage_error[j]; } PetscStackCallBLAS("LAPACKgetrs",LAPACKgetrs_("No transpose",&m,&one,ImV,&n,ipiv,scheme->gamma+1,&n,&info)); if (info < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Bad argument to GETRS"); if (info > 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Should not happen"); /* beta[0] (rho in B,J,W 2007) e.rho = 1/factorial(p+2) - glm.B(1,:)*glm.c.^(p+1)/factorial(p+1) ... + glm.V(1,2:end)*e.beta;% - e.epsilon; % Note: The paper (B,J,W 2007) includes the last term in their definition * */ scheme->beta[0] = 1./Factorial(p+2); for (j=0; jbeta[0] -= b[0*s+j]*CPowF(c[j],p+1); for (j=1; jbeta[0] += v[0*r+j]*scheme->beta[j]; /* gamma[0] (sigma in B,J,W 2007) * e.sigma = glm.B(1,:)*e.xi + glm.V(1,2:end)*e.gamma; * */ scheme->gamma[0] = 0.0; for (j=0; jgamma[0] += b[0*s+j]*scheme->stage_error[j]; for (j=1; jgamma[0] += v[0*s+j]*scheme->gamma[j]; /* Assemble H * % Determine the error estimators phi H = [[cpow(glm.c,p) + C*e.alpha] [cpow(glm.c,p+1) + C*e.beta] ... [e.xi - C*(e.gamma + 0*e.epsilon*eye(s-1,1))]]'; % Paper has formula above without the 0, but that term must be left % out to satisfy the conditions they propose and to make the % example schemes work e.H = H; e.phi = (H \ [1 0 0;1 1 0;0 0 -1])'; e.psi = -e.phi*C; * */ for (j=0; jstage_error[j]; for (k=1; kalpha[k]; H[1+j*3] += CPowF(c[j],k-1)*scheme->beta[k]; H[2+j*3] -= CPowF(c[j],k-1)*scheme->gamma[k]; } } bmat[0+0*ss] = 1.; bmat[0+1*ss] = 0.; bmat[0+2*ss] = 0.; bmat[1+0*ss] = 1.; bmat[1+1*ss] = 1.; bmat[1+2*ss] = 0.; bmat[2+0*ss] = 0.; bmat[2+1*ss] = 0.; bmat[2+2*ss] = -1.; m = 3; ierr = PetscBLASIntCast(s,&n);CHKERRQ(ierr); ierr = PetscBLASIntCast(ss,&ldb);CHKERRQ(ierr); rcond = 1e-12; #if defined(PETSC_MISSING_LAPACK_GELSS) /* ESSL does not have this routine */ SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GELSS - Lapack routine is unavailable\nNot able to run GL time stepping."); #else #if defined(PETSC_USE_COMPLEX) /* ZGELSS( M, N, NRHS, A, LDA, B, LDB, S, RCOND, RANK, WORK, LWORK, RWORK, INFO) */ PetscStackCallBLAS("LAPACKgelss",LAPACKgelss_(&m,&n,&m,H,&m,bmat,&ldb,sing,&rcond,&rank,workscalar,&lwork,workreal,&info)); #else /* DGELSS( M, N, NRHS, A, LDA, B, LDB, S, RCOND, RANK, WORK, LWORK, INFO) */ PetscStackCallBLAS("LAPACKgelss",LAPACKgelss_(&m,&n,&m,H,&m,bmat,&ldb,sing,&rcond,&rank,workscalar,&lwork,&info)); #endif if (info < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Bad argument to GELSS"); if (info > 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"SVD failed to converge"); #endif for (j=0; j<3; j++) { for (k=0; kphi[k+j*s] = bmat[k+j*ss]; } /* the other part of the error estimator, psi in B,J,W 2007 */ scheme->psi[0*r+0] = 0.; scheme->psi[1*r+0] = 0.; scheme->psi[2*r+0] = 0.; for (j=1; jpsi[0*r+j] = 0.; scheme->psi[1*r+j] = 0.; scheme->psi[2*r+j] = 0.; for (k=0; kpsi[0*r+j] -= CPowF(c[k],j-1)*scheme->phi[0*s+k]; scheme->psi[1*r+j] -= CPowF(c[k],j-1)*scheme->phi[1*s+k]; scheme->psi[2*r+j] -= CPowF(c[k],j-1)*scheme->phi[2*s+k]; } } ierr = PetscFree7(ImV,H,bmat,workscalar,workreal,sing,ipiv);CHKERRQ(ierr); } /* Check which properties are satisfied */ scheme->stiffly_accurate = PETSC_TRUE; if (scheme->c[s-1] != 1.) scheme->stiffly_accurate = PETSC_FALSE; for (j=0; jstiffly_accurate = PETSC_FALSE; for (j=0; jstiffly_accurate = PETSC_FALSE; scheme->fsal = scheme->stiffly_accurate; /* FSAL is stronger */ for (j=0; j1 && b[1*s+j] != 0.) scheme->fsal = PETSC_FALSE; if (b[1*s+r-1] != 1.) scheme->fsal = PETSC_FALSE; for (j=0; j1 && v[1*r+j] != 0.) scheme->fsal = PETSC_FALSE; *inscheme = scheme; PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "TSGLLESchemeDestroy" static PetscErrorCode TSGLLESchemeDestroy(TSGLLEScheme sc) { PetscErrorCode ierr; PetscFunctionBegin; ierr = PetscFree5(sc->c,sc->a,sc->b,sc->u,sc->v);CHKERRQ(ierr); ierr = PetscFree6(sc->alpha,sc->beta,sc->gamma,sc->phi,sc->psi,sc->stage_error);CHKERRQ(ierr); ierr = PetscFree(sc);CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "TSGLLEDestroy_Default" static PetscErrorCode TSGLLEDestroy_Default(TS_GLLE *gl) { PetscErrorCode ierr; PetscInt i; PetscFunctionBegin; for (i=0; inschemes; i++) { if (gl->schemes[i]) {ierr = TSGLLESchemeDestroy(gl->schemes[i]);CHKERRQ(ierr);} } ierr = PetscFree(gl->schemes);CHKERRQ(ierr); gl->nschemes = 0; ierr = PetscMemzero(gl->type_name,sizeof(gl->type_name));CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "TSGLLEViewTable_Private" static PetscErrorCode TSGLLEViewTable_Private(PetscViewer viewer,PetscInt m,PetscInt n,const PetscScalar a[],const char name[]) { PetscErrorCode ierr; PetscBool iascii; PetscInt i,j; PetscFunctionBegin; ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); if (iascii) { ierr = PetscViewerASCIIPrintf(viewer,"%30s = [",name);CHKERRQ(ierr); for (i=0; ip,sc->q,sc->r,sc->s);CHKERRQ(ierr); ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); ierr = PetscViewerASCIIPrintf(viewer,"Stiffly accurate: %s, FSAL: %s\n",sc->stiffly_accurate ? "yes" : "no",sc->fsal ? "yes" : "no");CHKERRQ(ierr); ierr = PetscViewerASCIIPrintf(viewer,"Leading error constants: %10.3e %10.3e %10.3e\n", PetscRealPart(sc->alpha[0]),PetscRealPart(sc->beta[0]),PetscRealPart(sc->gamma[0]));CHKERRQ(ierr); ierr = TSGLLEViewTable_Private(viewer,1,sc->s,sc->c,"Abscissas c");CHKERRQ(ierr); if (view_details) { ierr = TSGLLEViewTable_Private(viewer,sc->s,sc->s,sc->a,"A");CHKERRQ(ierr); ierr = TSGLLEViewTable_Private(viewer,sc->r,sc->s,sc->b,"B");CHKERRQ(ierr); ierr = TSGLLEViewTable_Private(viewer,sc->s,sc->r,sc->u,"U");CHKERRQ(ierr); ierr = TSGLLEViewTable_Private(viewer,sc->r,sc->r,sc->v,"V");CHKERRQ(ierr); ierr = TSGLLEViewTable_Private(viewer,3,sc->s,sc->phi,"Error estimate phi");CHKERRQ(ierr); ierr = TSGLLEViewTable_Private(viewer,3,sc->r,sc->psi,"Error estimate psi");CHKERRQ(ierr); ierr = TSGLLEViewTable_Private(viewer,1,sc->r,sc->alpha,"Modify alpha");CHKERRQ(ierr); ierr = TSGLLEViewTable_Private(viewer,1,sc->r,sc->beta,"Modify beta");CHKERRQ(ierr); ierr = TSGLLEViewTable_Private(viewer,1,sc->r,sc->gamma,"Modify gamma");CHKERRQ(ierr); ierr = TSGLLEViewTable_Private(viewer,1,sc->s,sc->stage_error,"Stage error xi");CHKERRQ(ierr); } ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Viewer type %s not supported",((PetscObject)viewer)->type_name); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "TSGLLEEstimateHigherMoments_Default" static PetscErrorCode TSGLLEEstimateHigherMoments_Default(TSGLLEScheme sc,PetscReal h,Vec Ydot[],Vec Xold[],Vec hm[]) { PetscErrorCode ierr; PetscInt i; PetscFunctionBegin; if (sc->r > 64 || sc->s > 64) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Ridiculous number of stages or items passed between stages"); /* build error vectors*/ for (i=0; i<3; i++) { PetscScalar phih[64]; PetscInt j; for (j=0; js; j++) phih[j] = sc->phi[i*sc->s+j]*h; ierr = VecZeroEntries(hm[i]);CHKERRQ(ierr); ierr = VecMAXPY(hm[i],sc->s,phih,Ydot);CHKERRQ(ierr); ierr = VecMAXPY(hm[i],sc->r,&sc->psi[i*sc->r],Xold);CHKERRQ(ierr); } PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "TSGLLECompleteStep_Rescale" static PetscErrorCode TSGLLECompleteStep_Rescale(TSGLLEScheme sc,PetscReal h,TSGLLEScheme next_sc,PetscReal next_h,Vec Ydot[],Vec Xold[],Vec X[]) { PetscErrorCode ierr; PetscScalar brow[32],vrow[32]; PetscInt i,j,r,s; PetscFunctionBegin; /* Build the new solution from (X,Ydot) */ r = sc->r; s = sc->s; for (i=0; ib[i*s+j]; ierr = VecMAXPY(X[i],s,brow,Ydot);CHKERRQ(ierr); for (j=0; jv[i*r+j]; ierr = VecMAXPY(X[i],r,vrow,Xold);CHKERRQ(ierr); } PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "TSGLLECompleteStep_RescaleAndModify" static PetscErrorCode TSGLLECompleteStep_RescaleAndModify(TSGLLEScheme sc,PetscReal h,TSGLLEScheme next_sc,PetscReal next_h,Vec Ydot[],Vec Xold[],Vec X[]) { PetscErrorCode ierr; PetscScalar brow[32],vrow[32]; PetscReal ratio; PetscInt i,j,p,r,s; PetscFunctionBegin; /* Build the new solution from (X,Ydot) */ p = sc->p; r = sc->r; s = sc->s; ratio = next_h/h; for (i=0; ib[i*s+j] + (PetscPowRealInt(ratio,i) - PetscPowRealInt(ratio,p+1))*(+ sc->alpha[i]*sc->phi[0*s+j]) + (PetscPowRealInt(ratio,i) - PetscPowRealInt(ratio,p+2))*(+ sc->beta [i]*sc->phi[1*s+j] + sc->gamma[i]*sc->phi[2*s+j])); } ierr = VecMAXPY(X[i],s,brow,Ydot);CHKERRQ(ierr); for (j=0; jv[i*r+j] + (PetscPowRealInt(ratio,i) - PetscPowRealInt(ratio,p+1))*(+ sc->alpha[i]*sc->psi[0*r+j]) + (PetscPowRealInt(ratio,i) - PetscPowRealInt(ratio,p+2))*(+ sc->beta [i]*sc->psi[1*r+j] + sc->gamma[i]*sc->psi[2*r+j])); } ierr = VecMAXPY(X[i],r,vrow,Xold);CHKERRQ(ierr); } if (r < next_sc->r) { if (r+1 != next_sc->r) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Cannot accommodate jump in r greater than 1"); ierr = VecZeroEntries(X[r]);CHKERRQ(ierr); for (j=0; jphi[0*s+j]; ierr = VecMAXPY(X[r],s,brow,Ydot);CHKERRQ(ierr); for (j=0; jpsi[0*r+j]; ierr = VecMAXPY(X[r],r,vrow,Xold);CHKERRQ(ierr); } PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "TSGLLECreate_IRKS" static PetscErrorCode TSGLLECreate_IRKS(TS ts) { TS_GLLE *gl = (TS_GLLE*)ts->data; PetscErrorCode ierr; PetscFunctionBegin; gl->Destroy = TSGLLEDestroy_Default; gl->EstimateHigherMoments = TSGLLEEstimateHigherMoments_Default; gl->CompleteStep = TSGLLECompleteStep_RescaleAndModify; ierr = PetscMalloc1(10,&gl->schemes);CHKERRQ(ierr); gl->nschemes = 0; { /* p=1,q=1, r=s=2, A- and L-stable with error estimates of order 2 and 3 * Listed in Butcher & Podhaisky 2006. On error estimation in general linear methods for stiff ODE. * irks(0.3,0,[.3,1],[1],1) * Note: can be made to have classical order (not stage order) 2 by replacing 0.3 with 1-sqrt(1/2) * but doing so would sacrifice the error estimator. */ const PetscScalar c[2] = {3./10., 1.}; const PetscScalar a[2][2] = {{3./10., 0}, {7./10., 3./10.}}; const PetscScalar b[2][2] = {{7./10., 3./10.}, {0,1}}; const PetscScalar u[2][2] = {{1,0},{1,0}}; const PetscScalar v[2][2] = {{1,0},{0,0}}; ierr = TSGLLESchemeCreate(1,1,2,2,c,*a,*b,*u,*v,&gl->schemes[gl->nschemes++]);CHKERRQ(ierr); } { /* p=q=2, r=s=3: irks(4/9,0,[1:3]/3,[0.33852],1) */ /* http://www.math.auckland.ac.nz/~hpod/atlas/i2a.html */ const PetscScalar c[3] = {1./3., 2./3., 1} ,a[3][3] = {{4./9. ,0 , 0}, {1.03750643704090e+00 , 4./9., 0}, {7.67024779410304e-01 , -3.81140216918943e-01, 4./9.}} ,b[3][3] = {{0.767024779410304, -0.381140216918943, 4./9.}, {0.000000000000000, 0.000000000000000, 1.000000000000000}, {-2.075048385225385, 0.621728385225383, 1.277197204924873}} ,u[3][3] = {{1.0000000000000000, -0.1111111111111109, -0.0925925925925922}, {1.0000000000000000, -0.8152842148186744, -0.4199095530877056}, {1.0000000000000000, 0.1696709930641948, 0.0539741070314165}} ,v[3][3] = {{1.0000000000000000, 0.1696709930641948, 0.0539741070314165}, {0.000000000000000, 0.000000000000000, 0.000000000000000}, {0.000000000000000, 0.176122795075129, 0.000000000000000}}; ierr = TSGLLESchemeCreate(2,2,3,3,c,*a,*b,*u,*v,&gl->schemes[gl->nschemes++]);CHKERRQ(ierr); } { /* p=q=3, r=s=4: irks(9/40,0,[1:4]/4,[0.3312 1.0050],[0.49541 1;1 0]) */ const PetscScalar c[4] = {0.25,0.5,0.75,1.0} ,a[4][4] = {{9./40. , 0, 0, 0}, {2.11286958887701e-01 , 9./40. , 0, 0}, {9.46338294287584e-01 , -3.42942861246094e-01, 9./40. , 0}, {0.521490453970721 , -0.662474225622980, 0.490476425459734, 9./40. }} ,b[4][4] = {{0.521490453970721 , -0.662474225622980, 0.490476425459734, 9./40. }, {0.000000000000000 , 0.000000000000000, 0.000000000000000, 1.000000000000000}, {-0.084677029310348 , 1.390757514776085, -1.568157386206001, 2.023192696767826}, {0.465383797936408 , 1.478273530625148, -1.930836081010182, 1.644872111193354}} ,u[4][4] = {{1.00000000000000000 , 0.02500000000001035, -0.02499999999999053, -0.00442708333332865}, {1.00000000000000000 , 0.06371304111232945, -0.04032173972189845, -0.01389438413189452}, {1.00000000000000000 , -0.07839543304147778, 0.04738685705116663, 0.02032603595928376}, {1.00000000000000000 , 0.42550734619251651, 0.10800718022400080, -0.01726712647760034}} ,v[4][4] = {{1.00000000000000000 , 0.42550734619251651, 0.10800718022400080, -0.01726712647760034}, {0.000000000000000 , 0.000000000000000, 0.000000000000000, 0.000000000000000}, {0.000000000000000 , -1.761115796027561, -0.521284157173780, 0.258249384305463}, {0.000000000000000 , -1.657693358744728, -1.052227765232394, 0.521284157173780}}; ierr = TSGLLESchemeCreate(3,3,4,4,c,*a,*b,*u,*v,&gl->schemes[gl->nschemes++]);CHKERRQ(ierr); } { /* p=q=4, r=s=5: irks(3/11,0,[1:5]/5, [0.1715 -0.1238 0.6617],... [ -0.0812 0.4079 1.0000 1.0000 0 0 0.8270 1.0000 0]) */ const PetscScalar c[5] = {0.2,0.4,0.6,0.8,1.0} ,a[5][5] = {{2.72727272727352e-01 , 0.00000000000000e+00, 0.00000000000000e+00 , 0.00000000000000e+00 , 0.00000000000000e+00}, {-1.03980153733431e-01, 2.72727272727405e-01, 0.00000000000000e+00, 0.00000000000000e+00 , 0.00000000000000e+00}, {-1.58615400341492e+00, 7.44168951881122e-01, 2.72727272727309e-01, 0.00000000000000e+00 , 0.00000000000000e+00}, {-8.73658042865628e-01, 5.37884671894595e-01, -1.63298538799523e-01, 2.72727272726996e-01 , 0.00000000000000e+00}, {2.95489397443992e-01 , -1.18481693910097e+00 , -6.68029812659953e-01 , 1.00716687860943e+00 , 2.72727272727288e-01}} ,b[5][5] = {{2.95489397443992e-01 , -1.18481693910097e+00 , -6.68029812659953e-01 , 1.00716687860943e+00 , 2.72727272727288e-01}, {0.00000000000000e+00 , 1.11022302462516e-16 , -2.22044604925031e-16 , 0.00000000000000e+00 , 1.00000000000000e+00}, {-4.05882503986005e+00, -4.00924006567769e+00, -1.38930610972481e+00, 4.45223930308488e+00 , 6.32331093108427e-01}, {8.35690179937017e+00 , -2.26640927349732e+00 , 6.86647884973826e+00 , -5.22595158025740e+00 , 4.50893068837431e+00}, {1.27656267027479e+01 , 2.80882153840821e+00 , 8.91173096522890e+00 , -1.07936444078906e+01 , 4.82534148988854e+00}} ,u[5][5] = {{1.00000000000000e+00 , -7.27272727273551e-02 , -3.45454545454419e-02 , -4.12121212119565e-03 ,-2.96969696964014e-04}, {1.00000000000000e+00 , 2.31252881006154e-01 , -8.29487834416481e-03 , -9.07191207681020e-03 ,-1.70378403743473e-03}, {1.00000000000000e+00 , 1.16925777880663e+00 , 3.59268562942635e-02 , -4.09013451730615e-02 ,-1.02411119670164e-02}, {1.00000000000000e+00 , 1.02634463704356e+00 , 1.59375044913405e-01 , 1.89673015035370e-03 ,-4.89987231897569e-03}, {1.00000000000000e+00 , 1.27746320298021e+00 , 2.37186008132728e-01 , -8.28694373940065e-02 ,-5.34396510196430e-02}} ,v[5][5] = {{1.00000000000000e+00 , 1.27746320298021e+00 , 2.37186008132728e-01 , -8.28694373940065e-02 ,-5.34396510196430e-02}, {0.00000000000000e+00 , -1.77635683940025e-15 , -1.99840144432528e-15 , -9.99200722162641e-16 ,-3.33066907387547e-16}, {0.00000000000000e+00 , 4.37280081906924e+00 , 5.49221645016377e-02 , -8.88913177394943e-02 , 1.12879077989154e-01}, {0.00000000000000e+00 , -1.22399504837280e+01 , -5.21287338448645e+00 , -8.03952325565291e-01 , 4.60298678047147e-01}, {0.00000000000000e+00 , -1.85178762883829e+01 , -5.21411849862624e+00 , -1.04283436528809e+00 , 7.49030161063651e-01}}; ierr = TSGLLESchemeCreate(4,4,5,5,c,*a,*b,*u,*v,&gl->schemes[gl->nschemes++]);CHKERRQ(ierr); } { /* p=q=5, r=s=6; irks(1/3,0,[1:6]/6,... [-0.0489 0.4228 -0.8814 0.9021],... [-0.3474 -0.6617 0.6294 0.2129 0.0044 -0.4256 -0.1427 -0.8936 -0.8267 0.4821 0.1371 -0.2557 -0.4426 -0.3855 -0.7514 0.3014]) */ const PetscScalar c[6] = {1./6, 2./6, 3./6, 4./6, 5./6, 1.} ,a[6][6] = {{ 3.33333333333940e-01, 0 , 0 , 0 , 0 , 0 }, { -8.64423857333350e-02, 3.33333333332888e-01, 0 , 0 , 0 , 0 }, { -2.16850174258252e+00, -2.23619072028839e+00, 3.33333333335204e-01, 0 , 0 , 0 }, { -4.73160970138997e+00, -3.89265344629268e+00, -2.76318716520933e-01, 3.33333333335759e-01, 0 , 0 }, { -6.75187540297338e+00, -7.90756533769377e+00, 7.90245051802259e-01, -4.48352364517632e-01, 3.33333333328483e-01, 0 }, { -4.26488287921548e+00, -1.19320395589302e+01, 3.38924509887755e+00, -2.23969848002481e+00, 6.62807710124007e-01, 3.33333333335440e-01}} ,b[6][6] = {{ -4.26488287921548e+00, -1.19320395589302e+01, 3.38924509887755e+00, -2.23969848002481e+00, 6.62807710124007e-01, 3.33333333335440e-01}, { -8.88178419700125e-16, 4.44089209850063e-16, -1.54737334057131e-15, -8.88178419700125e-16, 0.00000000000000e+00, 1.00000000000001e+00}, { -2.87780425770651e+01, -1.13520448264971e+01, 2.62002318943161e+01, 2.56943874812797e+01, -3.06702268304488e+01, 6.68067773510103e+00}, { 5.47971245256474e+01, 6.80366875868284e+01, -6.50952588861999e+01, -8.28643975339097e+01, 8.17416943896414e+01, -1.17819043489036e+01}, { -2.33332114788869e+02, 6.12942539462634e+01, -4.91850135865944e+01, 1.82716844135480e+02, -1.29788173979395e+02, 3.09968095651099e+01}, { -1.72049132343751e+02, 8.60194713593999e+00, 7.98154219170200e-01, 1.50371386053218e+02, -1.18515423962066e+02, 2.50898277784663e+01}} ,u[6][6] = {{ 1.00000000000000e+00, -1.66666666666870e-01, -4.16666666664335e-02, -3.85802469124815e-03, -2.25051440302250e-04, -9.64506172339142e-06}, { 1.00000000000000e+00, 8.64423857327162e-02, -4.11484912671353e-02, -1.11450903217645e-02, -1.47651050487126e-03, -1.34395070766826e-04}, { 1.00000000000000e+00, 4.57135912953434e+00, 1.06514719719137e+00, 1.33517564218007e-01, 1.11365952968659e-02, 6.12382756769504e-04}, { 1.00000000000000e+00, 9.23391519753404e+00, 2.22431212392095e+00, 2.91823807741891e-01, 2.52058456411084e-02, 1.22800542949647e-03}, { 1.00000000000000e+00, 1.48175480533865e+01, 3.73439117461835e+00, 5.14648336541804e-01, 4.76430038853402e-02, 2.56798515502156e-03}, { 1.00000000000000e+00, 1.50512347758335e+01, 4.10099701165164e+00, 5.66039141003603e-01, 3.91213893800891e-02, -2.99136269067853e-03}} ,v[6][6] = {{ 1.00000000000000e+00, 1.50512347758335e+01, 4.10099701165164e+00, 5.66039141003603e-01, 3.91213893800891e-02, -2.99136269067853e-03}, { 0.00000000000000e+00, -4.88498130835069e-15, -6.43929354282591e-15, -3.55271367880050e-15, -1.22124532708767e-15, -3.12250225675825e-16}, { 0.00000000000000e+00, 1.22250171233141e+01, -1.77150760606169e+00, 3.54516769879390e-01, 6.22298845883398e-01, 2.31647447450276e-01}, { 0.00000000000000e+00, -4.48339457331040e+01, -3.57363126641880e-01, 5.18750173123425e-01, 6.55727990241799e-02, 1.63175368287079e-01}, { 0.00000000000000e+00, 1.37297394708005e+02, -1.60145272991317e+00, -5.05319555199441e+00, 1.55328940390990e-01, 9.16629423682464e-01}, { 0.00000000000000e+00, 1.05703241119022e+02, -1.16610260983038e+00, -2.99767252773859e+00, -1.13472315553890e-01, 1.09742849254729e+00}}; ierr = TSGLLESchemeCreate(5,5,6,6,c,*a,*b,*u,*v,&gl->schemes[gl->nschemes++]);CHKERRQ(ierr); } PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "TSGLLESetType" /*@C TSGLLESetType - sets the class of general linear method to use for time-stepping Collective on TS Input Parameters: + ts - the TS context - type - a method Options Database Key: . -ts_gl_type - sets the method, use -help for a list of available method (e.g. irks) Notes: See "petsc/include/petscts.h" for available methods (for instance) . TSGLLE_IRKS - Diagonally implicit methods with inherent Runge-Kutta stability (for stiff problems) Normally, it is best to use the TSSetFromOptions() command and then set the TSGLLE type from the options database rather than by using this routine. Using the options database provides the user with maximum flexibility in evaluating the many different solvers. The TSGLLESetType() routine is provided for those situations where it is necessary to set the timestepping solver independently of the command line or options database. This might be the case, for example, when the choice of solver changes during the execution of the program, and the user's application is taking responsibility for choosing the appropriate method. Level: intermediate .keywords: TS, TSGLLE, set, type @*/ PetscErrorCode TSGLLESetType(TS ts,TSGLLEType type) { PetscErrorCode ierr; PetscFunctionBegin; PetscValidHeaderSpecific(ts,TS_CLASSID,1); PetscValidCharPointer(type,2); ierr = PetscTryMethod(ts,"TSGLLESetType_C",(TS,TSGLLEType),(ts,type));CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "TSGLLESetAcceptType" /*@C TSGLLESetAcceptType - sets the acceptance test Time integrators that need to control error must have the option to reject a time step based on local error estimates. This function allows different schemes to be set. Logically Collective on TS Input Parameters: + ts - the TS context - type - the type Options Database Key: . -ts_gl_accept_type - sets the method used to determine whether to accept or reject a step Level: intermediate .seealso: TS, TSGLLE, TSGLLEAcceptRegister(), TSGLLEAdapt, set type @*/ PetscErrorCode TSGLLESetAcceptType(TS ts,TSGLLEAcceptType type) { PetscErrorCode ierr; PetscFunctionBegin; PetscValidHeaderSpecific(ts,TS_CLASSID,1); PetscValidCharPointer(type,2); ierr = PetscTryMethod(ts,"TSGLLESetAcceptType_C",(TS,TSGLLEAcceptType),(ts,type));CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "TSGLLEGetAdapt" /*@C TSGLLEGetAdapt - gets the TSGLLEAdapt object from the TS Not Collective Input Parameter: . ts - the TS context Output Parameter: . adapt - the TSGLLEAdapt context Notes: This allows the user set options on the TSGLLEAdapt object. Usually it is better to do this using the options database, so this function is rarely needed. Level: advanced .seealso: TSGLLEAdapt, TSGLLEAdaptRegister() @*/ PetscErrorCode TSGLLEGetAdapt(TS ts,TSGLLEAdapt *adapt) { PetscErrorCode ierr; PetscFunctionBegin; PetscValidHeaderSpecific(ts,TS_CLASSID,1); PetscValidPointer(adapt,2); ierr = PetscUseMethod(ts,"TSGLLEGetAdapt_C",(TS,TSGLLEAdapt*),(ts,adapt));CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "TSGLLEAccept_Always" static PetscErrorCode TSGLLEAccept_Always(TS ts,PetscReal tleft,PetscReal h,const PetscReal enorms[],PetscBool *accept) { PetscFunctionBegin; *accept = PETSC_TRUE; PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "TSGLLEUpdateWRMS" static PetscErrorCode TSGLLEUpdateWRMS(TS ts) { TS_GLLE *gl = (TS_GLLE*)ts->data; PetscErrorCode ierr; PetscScalar *x,*w; PetscInt n,i; PetscFunctionBegin; ierr = VecGetArray(gl->X[0],&x);CHKERRQ(ierr); ierr = VecGetArray(gl->W,&w);CHKERRQ(ierr); ierr = VecGetLocalSize(gl->W,&n);CHKERRQ(ierr); for (i=0; iwrms_atol + gl->wrms_rtol*PetscAbsScalar(x[i])); ierr = VecRestoreArray(gl->X[0],&x);CHKERRQ(ierr); ierr = VecRestoreArray(gl->W,&w);CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "TSGLLEVecNormWRMS" static PetscErrorCode TSGLLEVecNormWRMS(TS ts,Vec X,PetscReal *nrm) { TS_GLLE *gl = (TS_GLLE*)ts->data; PetscErrorCode ierr; PetscScalar *x,*w; PetscReal sum = 0.0,gsum; PetscInt n,N,i; PetscFunctionBegin; ierr = VecGetArray(X,&x);CHKERRQ(ierr); ierr = VecGetArray(gl->W,&w);CHKERRQ(ierr); ierr = VecGetLocalSize(gl->W,&n);CHKERRQ(ierr); for (i=0; iW,&w);CHKERRQ(ierr); ierr = MPIU_Allreduce(&sum,&gsum,1,MPIU_REAL,MPIU_SUM,PetscObjectComm((PetscObject)ts));CHKERRQ(ierr); ierr = VecGetSize(gl->W,&N);CHKERRQ(ierr); *nrm = PetscSqrtReal(gsum/(1.*N)); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "TSGLLESetType_GLLE" static PetscErrorCode TSGLLESetType_GLLE(TS ts,TSGLLEType type) { PetscErrorCode ierr,(*r)(TS); PetscBool same; TS_GLLE *gl = (TS_GLLE*)ts->data; PetscFunctionBegin; if (gl->type_name[0]) { ierr = PetscStrcmp(gl->type_name,type,&same);CHKERRQ(ierr); if (same) PetscFunctionReturn(0); ierr = (*gl->Destroy)(gl);CHKERRQ(ierr); } ierr = PetscFunctionListFind(TSGLLEList,type,&r);CHKERRQ(ierr); if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown TSGLLE type \"%s\" given",type); ierr = (*r)(ts);CHKERRQ(ierr); ierr = PetscStrcpy(gl->type_name,type);CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "TSGLLESetAcceptType_GLLE" static PetscErrorCode TSGLLESetAcceptType_GLLE(TS ts,TSGLLEAcceptType type) { PetscErrorCode ierr; TSGLLEAcceptFunction r; TS_GLLE *gl = (TS_GLLE*)ts->data; PetscFunctionBegin; ierr = PetscFunctionListFind(TSGLLEAcceptList,type,&r);CHKERRQ(ierr); if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown TSGLLEAccept type \"%s\" given",type); gl->Accept = r; ierr = PetscStrncpy(gl->accept_name,type,sizeof(gl->accept_name));CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "TSGLLEGetAdapt_GLLE" static PetscErrorCode TSGLLEGetAdapt_GLLE(TS ts,TSGLLEAdapt *adapt) { PetscErrorCode ierr; TS_GLLE *gl = (TS_GLLE*)ts->data; PetscFunctionBegin; if (!gl->adapt) { ierr = TSGLLEAdaptCreate(PetscObjectComm((PetscObject)ts),&gl->adapt);CHKERRQ(ierr); ierr = PetscObjectIncrementTabLevel((PetscObject)gl->adapt,(PetscObject)ts,1);CHKERRQ(ierr); ierr = PetscLogObjectParent((PetscObject)ts,(PetscObject)gl->adapt);CHKERRQ(ierr); } *adapt = gl->adapt; PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "TSGLLEChooseNextScheme" static PetscErrorCode TSGLLEChooseNextScheme(TS ts,PetscReal h,const PetscReal hmnorm[],PetscInt *next_scheme,PetscReal *next_h,PetscBool *finish) { PetscErrorCode ierr; TS_GLLE *gl = (TS_GLLE*)ts->data; PetscInt i,n,cur_p,cur,next_sc,candidates[64],orders[64]; PetscReal errors[64],costs[64],tleft; PetscFunctionBegin; cur = -1; cur_p = gl->schemes[gl->current_scheme]->p; tleft = ts->max_time - (ts->ptime + ts->time_step); for (i=0,n=0; inschemes; i++) { TSGLLEScheme sc = gl->schemes[i]; if (sc->p < gl->min_order || gl->max_order < sc->p) continue; if (sc->p == cur_p - 1) errors[n] = PetscAbsScalar(sc->alpha[0])*hmnorm[0]; else if (sc->p == cur_p) errors[n] = PetscAbsScalar(sc->alpha[0])*hmnorm[1]; else if (sc->p == cur_p+1) errors[n] = PetscAbsScalar(sc->alpha[0])*(hmnorm[2]+hmnorm[3]); else continue; candidates[n] = i; orders[n] = PetscMin(sc->p,sc->q); /* order of global truncation error */ costs[n] = sc->s; /* estimate the cost as the number of stages */ if (i == gl->current_scheme) cur = n; n++; } if (cur < 0 || gl->nschemes <= cur) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Current scheme not found in scheme list"); ierr = TSGLLEAdaptChoose(gl->adapt,n,orders,errors,costs,cur,h,tleft,&next_sc,next_h,finish);CHKERRQ(ierr); *next_scheme = candidates[next_sc]; ierr = PetscInfo7(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);CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "TSGLLEGetMaxSizes" static PetscErrorCode TSGLLEGetMaxSizes(TS ts,PetscInt *max_r,PetscInt *max_s) { TS_GLLE *gl = (TS_GLLE*)ts->data; PetscFunctionBegin; *max_r = gl->schemes[gl->nschemes-1]->r; *max_s = gl->schemes[gl->nschemes-1]->s; PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "TSSolve_GLLE" static PetscErrorCode TSSolve_GLLE(TS ts) { TS_GLLE *gl = (TS_GLLE*)ts->data; PetscInt i,k,its,lits,max_r,max_s; PetscBool final_step,finish; SNESConvergedReason snesreason; PetscErrorCode ierr; PetscFunctionBegin; ierr = TSMonitor(ts,ts->steps,ts->ptime,ts->vec_sol);CHKERRQ(ierr); ierr = TSGLLEGetMaxSizes(ts,&max_r,&max_s);CHKERRQ(ierr); ierr = VecCopy(ts->vec_sol,gl->X[0]);CHKERRQ(ierr); for (i=1; iX[i]);CHKERRQ(ierr); } ierr = TSGLLEUpdateWRMS(ts);CHKERRQ(ierr); if (0) { /* Find consistent initial data for DAE */ gl->stage_time = ts->ptime + ts->time_step; gl->scoeff = 1.; gl->stage = 0; ierr = VecCopy(ts->vec_sol,gl->Z);CHKERRQ(ierr); ierr = VecCopy(ts->vec_sol,gl->Y);CHKERRQ(ierr); ierr = SNESSolve(ts->snes,NULL,gl->Y);CHKERRQ(ierr); ierr = SNESGetIterationNumber(ts->snes,&its);CHKERRQ(ierr); ierr = SNESGetLinearSolveIterations(ts->snes,&lits);CHKERRQ(ierr); ierr = SNESGetConvergedReason(ts->snes,&snesreason);CHKERRQ(ierr); ts->snes_its += its; ts->ksp_its += lits; if (snesreason < 0 && ts->max_snes_failures > 0 && ++ts->num_snes_failures >= ts->max_snes_failures) { ts->reason = TS_DIVERGED_NONLINEAR_SOLVE; ierr = PetscInfo2(ts,"Step=%D, nonlinear solve solve failures %D greater than current TS allowed, stopping solve\n",ts->steps,ts->num_snes_failures);CHKERRQ(ierr); PetscFunctionReturn(0); } } if (gl->current_scheme < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"A starting scheme has not been provided"); for (k=0,final_step=PETSC_FALSE,finish=PETSC_FALSE; kmax_steps && !finish; k++) { PetscInt j,r,s,next_scheme = 0,rejections; PetscReal h,hmnorm[4],enorm[3],next_h; PetscBool accept; const PetscScalar *c,*a,*u; Vec *X,*Ydot,Y; TSGLLEScheme scheme = gl->schemes[gl->current_scheme]; r = scheme->r; s = scheme->s; c = scheme->c; a = scheme->a; u = scheme->u; h = ts->time_step; X = gl->X; Ydot = gl->Ydot; Y = gl->Y; if (ts->ptime > ts->max_time) break; /* We only call PreStep at the start of each STEP, not each STAGE. This is because it is possible to fail (have to restart a step) after multiple stages. */ ierr = TSPreStep(ts);CHKERRQ(ierr); rejections = 0; while (1) { for (i=0; iscoeff = 1./PetscRealPart(a[i*s+i]); shift = gl->scoeff/ts->time_step; gl->stage = i; gl->stage_time = ts->ptime + PetscRealPart(c[i])*h; /* * Stage equation: Y = h A Y' + U X * We assume that A is lower-triangular so that we can solve the stages (Y,Y') sequentially * Build the affine vector z_i = -[1/(h a_ii)](h sum_j a_ij y'_j + sum_j u_ij x_j) * Then y'_i = z + 1/(h a_ii) y_i */ ierr = VecZeroEntries(gl->Z);CHKERRQ(ierr); for (j=0; jZ,-shift*u[i*r+j],X[j]);CHKERRQ(ierr); } for (j=0; jZ,-shift*h*a[i*s+j],Ydot[j]);CHKERRQ(ierr); } /* Note: Z is used within function evaluation, Ydot = Z + shift*Y */ /* Compute an estimate of Y to start Newton iteration */ if (gl->extrapolate) { if (i==0) { /* Linear extrapolation on the first stage */ ierr = VecWAXPY(Y,c[i]*h,X[1],X[0]);CHKERRQ(ierr); } else { /* Linear extrapolation from the last stage */ ierr = VecAXPY(Y,(c[i]-c[i-1])*h,Ydot[i-1]);CHKERRQ(ierr); } } else if (i==0) { /* Directly use solution from the last step, otherwise reuse the last stage (do nothing) */ ierr = VecCopy(X[0],Y);CHKERRQ(ierr); } /* Solve this stage (Ydot[i] is computed during function evaluation) */ ierr = SNESSolve(ts->snes,NULL,Y);CHKERRQ(ierr); ierr = SNESGetIterationNumber(ts->snes,&its);CHKERRQ(ierr); ierr = SNESGetLinearSolveIterations(ts->snes,&lits);CHKERRQ(ierr); ierr = SNESGetConvergedReason(ts->snes,&snesreason);CHKERRQ(ierr); ts->snes_its += its; ts->ksp_its += lits; if (snesreason < 0 && ts->max_snes_failures > 0 && ++ts->num_snes_failures >= ts->max_snes_failures) { ts->reason = TS_DIVERGED_NONLINEAR_SOLVE; ierr = PetscInfo2(ts,"Step=%D, nonlinear solve solve failures %D greater than current TS allowed, stopping solve\n",ts->steps,ts->num_snes_failures);CHKERRQ(ierr); PetscFunctionReturn(0); } } gl->stage_time = ts->ptime + ts->time_step; ierr = (*gl->EstimateHigherMoments)(scheme,h,Ydot,gl->X,gl->himom);CHKERRQ(ierr); /* hmnorm[i] = h^{p+i}x^{(p+i)} with i=0,1,2; hmnorm[3] = h^{p+2}(dx'/dx) x^{(p+1)} */ for (i=0; i<3; i++) { ierr = TSGLLEVecNormWRMS(ts,gl->himom[i],&hmnorm[i+1]);CHKERRQ(ierr); } enorm[0] = PetscRealPart(scheme->alpha[0])*hmnorm[1]; enorm[1] = PetscRealPart(scheme->beta[0]) *hmnorm[2]; enorm[2] = PetscRealPart(scheme->gamma[0])*hmnorm[3]; ierr = (*gl->Accept)(ts,ts->max_time-gl->stage_time,h,enorm,&accept);CHKERRQ(ierr); if (accept) goto accepted; rejections++; ierr = PetscInfo3(ts,"Step %D (t=%g) not accepted, rejections=%D\n",k,gl->stage_time,rejections);CHKERRQ(ierr); if (rejections > gl->max_step_rejections) break; /* There are lots of reasons why a step might be rejected, including solvers not converging and other factors that TSGLLEChooseNextScheme does not support. Additionally, the error estimates may be very screwed up, so I'm not convinced that it's safe to just compute a new error estimate using the same interface as the current adaptor (the adaptor interface probably has to change). Here we make an arbitrary and naive choice. This assumes that steps were written in Nordsieck form. The "correct" method would be to re-complete the previous time step with the correct "next" step size. It is unclear to me whether the present ad-hoc method of rescaling X is stable. */ h *= 0.5; for (i=1; ir; i++) { ierr = VecScale(X[i],PetscPowRealInt(0.5,i));CHKERRQ(ierr); } } SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_CONV_FAILED,"Time step %D (t=%g) not accepted after %D failures\n",k,gl->stage_time,rejections); accepted: /* This term is not error, but it *would* be the leading term for a lower order method */ ierr = TSGLLEVecNormWRMS(ts,gl->X[scheme->r-1],&hmnorm[0]);CHKERRQ(ierr); /* Correct scaling so that these are equivalent to norms of the Nordsieck vectors */ ierr = PetscInfo4(ts,"Last moment norm %10.2e, estimated error norms %10.2e %10.2e %10.2e\n",hmnorm[0],enorm[0],enorm[1],enorm[2]);CHKERRQ(ierr); if (!final_step) { ierr = TSGLLEChooseNextScheme(ts,h,hmnorm,&next_scheme,&next_h,&final_step);CHKERRQ(ierr); } else { /* Dummy values to complete the current step in a consistent manner */ next_scheme = gl->current_scheme; next_h = h; finish = PETSC_TRUE; } X = gl->Xold; gl->Xold = gl->X; gl->X = X; ierr = (*gl->CompleteStep)(scheme,h,gl->schemes[next_scheme],next_h,Ydot,gl->Xold,gl->X);CHKERRQ(ierr); ierr = TSGLLEUpdateWRMS(ts);CHKERRQ(ierr); /* Post the solution for the user, we could avoid this copy with a small bit of cleverness */ ierr = VecCopy(gl->X[0],ts->vec_sol);CHKERRQ(ierr); ts->ptime += h; ts->steps++; ts->total_steps++; ierr = TSPostEvaluate(ts);CHKERRQ(ierr); ierr = TSPostStep(ts);CHKERRQ(ierr); ierr = TSMonitor(ts,ts->steps,ts->ptime,ts->vec_sol);CHKERRQ(ierr); gl->current_scheme = next_scheme; ts->time_step = next_h; } PetscFunctionReturn(0); } /*------------------------------------------------------------*/ #undef __FUNCT__ #define __FUNCT__ "TSReset_GLLE" static PetscErrorCode TSReset_GLLE(TS ts) { TS_GLLE *gl = (TS_GLLE*)ts->data; PetscInt max_r,max_s; PetscErrorCode ierr; PetscFunctionBegin; if (gl->setupcalled) { ierr = TSGLLEGetMaxSizes(ts,&max_r,&max_s);CHKERRQ(ierr); ierr = VecDestroyVecs(max_r,&gl->Xold);CHKERRQ(ierr); ierr = VecDestroyVecs(max_r,&gl->X);CHKERRQ(ierr); ierr = VecDestroyVecs(max_s,&gl->Ydot);CHKERRQ(ierr); ierr = VecDestroyVecs(3,&gl->himom);CHKERRQ(ierr); ierr = VecDestroy(&gl->W);CHKERRQ(ierr); ierr = VecDestroy(&gl->Y);CHKERRQ(ierr); ierr = VecDestroy(&gl->Z);CHKERRQ(ierr); } gl->setupcalled = PETSC_FALSE; PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "TSDestroy_GLLE" static PetscErrorCode TSDestroy_GLLE(TS ts) { TS_GLLE *gl = (TS_GLLE*)ts->data; PetscErrorCode ierr; PetscFunctionBegin; ierr = TSReset_GLLE(ts);CHKERRQ(ierr); if (gl->adapt) {ierr = TSGLLEAdaptDestroy(&gl->adapt);CHKERRQ(ierr);} if (gl->Destroy) {ierr = (*gl->Destroy)(gl);CHKERRQ(ierr);} ierr = PetscFree(ts->data);CHKERRQ(ierr); ierr = PetscObjectComposeFunction((PetscObject)ts,"TSGLLESetType_C",NULL);CHKERRQ(ierr); ierr = PetscObjectComposeFunction((PetscObject)ts,"TSGLLESetAcceptType_C",NULL);CHKERRQ(ierr); ierr = PetscObjectComposeFunction((PetscObject)ts,"TSGLLEGetAdapt_C",NULL);CHKERRQ(ierr); PetscFunctionReturn(0); } /* This defines the nonlinear equation that is to be solved with SNES g(x) = f(t,x,z+shift*x) = 0 */ #undef __FUNCT__ #define __FUNCT__ "SNESTSFormFunction_GLLE" static PetscErrorCode SNESTSFormFunction_GLLE(SNES snes,Vec x,Vec f,TS ts) { TS_GLLE *gl = (TS_GLLE*)ts->data; PetscErrorCode ierr; Vec Z,Ydot; DM dm,dmsave; PetscFunctionBegin; ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); ierr = TSGLLEGetVecs(ts,dm,&Z,&Ydot);CHKERRQ(ierr); ierr = VecWAXPY(Ydot,gl->scoeff/ts->time_step,x,Z);CHKERRQ(ierr); dmsave = ts->dm; ts->dm = dm; ierr = TSComputeIFunction(ts,gl->stage_time,x,Ydot,f,PETSC_FALSE);CHKERRQ(ierr); ts->dm = dmsave; ierr = TSGLLERestoreVecs(ts,dm,&Z,&Ydot);CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "SNESTSFormJacobian_GLLE" static PetscErrorCode SNESTSFormJacobian_GLLE(SNES snes,Vec x,Mat A,Mat B,TS ts) { TS_GLLE *gl = (TS_GLLE*)ts->data; PetscErrorCode ierr; Vec Z,Ydot; DM dm,dmsave; PetscFunctionBegin; ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); ierr = TSGLLEGetVecs(ts,dm,&Z,&Ydot);CHKERRQ(ierr); dmsave = ts->dm; ts->dm = dm; /* gl->Xdot will have already been computed in SNESTSFormFunction_GLLE */ ierr = TSComputeIJacobian(ts,gl->stage_time,x,gl->Ydot[gl->stage],gl->scoeff/ts->time_step,A,B,PETSC_FALSE);CHKERRQ(ierr); ts->dm = dmsave; ierr = TSGLLERestoreVecs(ts,dm,&Z,&Ydot);CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "TSSetUp_GLLE" static PetscErrorCode TSSetUp_GLLE(TS ts) { TS_GLLE *gl = (TS_GLLE*)ts->data; PetscInt max_r,max_s; PetscErrorCode ierr; DM dm; PetscFunctionBegin; gl->setupcalled = PETSC_TRUE; ierr = TSGLLEGetMaxSizes(ts,&max_r,&max_s);CHKERRQ(ierr); ierr = VecDuplicateVecs(ts->vec_sol,max_r,&gl->X);CHKERRQ(ierr); ierr = VecDuplicateVecs(ts->vec_sol,max_r,&gl->Xold);CHKERRQ(ierr); ierr = VecDuplicateVecs(ts->vec_sol,max_s,&gl->Ydot);CHKERRQ(ierr); ierr = VecDuplicateVecs(ts->vec_sol,3,&gl->himom);CHKERRQ(ierr); ierr = VecDuplicate(ts->vec_sol,&gl->W);CHKERRQ(ierr); ierr = VecDuplicate(ts->vec_sol,&gl->Y);CHKERRQ(ierr); ierr = VecDuplicate(ts->vec_sol,&gl->Z);CHKERRQ(ierr); /* Default acceptance tests and adaptivity */ if (!gl->Accept) {ierr = TSGLLESetAcceptType(ts,TSGLLEACCEPT_ALWAYS);CHKERRQ(ierr);} if (!gl->adapt) {ierr = TSGLLEGetAdapt(ts,&gl->adapt);CHKERRQ(ierr);} if (gl->current_scheme < 0) { PetscInt i; for (i=0;; i++) { if (gl->schemes[i]->p == gl->start_order) break; if (i+1 == gl->nschemes) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"No schemes available with requested start order %d",i); } gl->current_scheme = i; } ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); if (dm) { ierr = DMCoarsenHookAdd(dm,DMCoarsenHook_TSGLLE,DMRestrictHook_TSGLLE,ts);CHKERRQ(ierr); ierr = DMSubDomainHookAdd(dm,DMSubDomainHook_TSGLLE,DMSubDomainRestrictHook_TSGLLE,ts);CHKERRQ(ierr); } PetscFunctionReturn(0); } /*------------------------------------------------------------*/ #undef __FUNCT__ #define __FUNCT__ "TSSetFromOptions_GLLE" static PetscErrorCode TSSetFromOptions_GLLE(PetscOptionItems *PetscOptionsObject,TS ts) { TS_GLLE *gl = (TS_GLLE*)ts->data; char tname[256] = TSGLLE_IRKS,completef[256] = "rescale-and-modify"; PetscErrorCode ierr; PetscFunctionBegin; ierr = PetscOptionsHead(PetscOptionsObject,"General Linear ODE solver options");CHKERRQ(ierr); { PetscBool flg; ierr = PetscOptionsFList("-ts_gl_type","Type of GL method","TSGLLESetType",TSGLLEList,gl->type_name[0] ? gl->type_name : tname,tname,sizeof(tname),&flg);CHKERRQ(ierr); if (flg || !gl->type_name[0]) { ierr = TSGLLESetType(ts,tname);CHKERRQ(ierr); } ierr = PetscOptionsInt("-ts_gl_max_step_rejections","Maximum number of times to attempt a step","None",gl->max_step_rejections,&gl->max_step_rejections,NULL);CHKERRQ(ierr); ierr = PetscOptionsInt("-ts_gl_max_order","Maximum order to try","TSGLLESetMaxOrder",gl->max_order,&gl->max_order,NULL);CHKERRQ(ierr); ierr = PetscOptionsInt("-ts_gl_min_order","Minimum order to try","TSGLLESetMinOrder",gl->min_order,&gl->min_order,NULL);CHKERRQ(ierr); ierr = PetscOptionsInt("-ts_gl_start_order","Initial order to try","TSGLLESetMinOrder",gl->start_order,&gl->start_order,NULL);CHKERRQ(ierr); ierr = PetscOptionsEnum("-ts_gl_error_direction","Which direction to look when estimating error","TSGLLESetErrorDirection",TSGLLEErrorDirections,(PetscEnum)gl->error_direction,(PetscEnum*)&gl->error_direction,NULL);CHKERRQ(ierr); ierr = PetscOptionsBool("-ts_gl_extrapolate","Extrapolate stage solution from previous solution (sometimes unstable)","TSGLLESetExtrapolate",gl->extrapolate,&gl->extrapolate,NULL);CHKERRQ(ierr); ierr = PetscOptionsReal("-ts_gl_atol","Absolute tolerance","TSGLLESetTolerances",gl->wrms_atol,&gl->wrms_atol,NULL);CHKERRQ(ierr); ierr = PetscOptionsReal("-ts_gl_rtol","Relative tolerance","TSGLLESetTolerances",gl->wrms_rtol,&gl->wrms_rtol,NULL);CHKERRQ(ierr); ierr = PetscOptionsString("-ts_gl_complete","Method to use for completing the step","none",completef,completef,sizeof(completef),&flg);CHKERRQ(ierr); if (flg) { PetscBool match1,match2; ierr = PetscStrcmp(completef,"rescale",&match1);CHKERRQ(ierr); ierr = PetscStrcmp(completef,"rescale-and-modify",&match2);CHKERRQ(ierr); if (match1) gl->CompleteStep = TSGLLECompleteStep_Rescale; else if (match2) gl->CompleteStep = TSGLLECompleteStep_RescaleAndModify; else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"%s",completef); } { char type[256] = TSGLLEACCEPT_ALWAYS; ierr = 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);CHKERRQ(ierr); if (flg || !gl->accept_name[0]) { ierr = TSGLLESetAcceptType(ts,type);CHKERRQ(ierr); } } { TSGLLEAdapt adapt; ierr = TSGLLEGetAdapt(ts,&adapt);CHKERRQ(ierr); ierr = TSGLLEAdaptSetFromOptions(PetscOptionsObject,adapt);CHKERRQ(ierr); } } ierr = PetscOptionsTail();CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "TSView_GLLE" static PetscErrorCode TSView_GLLE(TS ts,PetscViewer viewer) { TS_GLLE *gl = (TS_GLLE*)ts->data; PetscInt i; PetscBool iascii,details; PetscErrorCode ierr; PetscFunctionBegin; ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); if (iascii) { ierr = PetscViewerASCIIPrintf(viewer," min order %D, max order %D, current order %D\n",gl->min_order,gl->max_order,gl->schemes[gl->current_scheme]->p);CHKERRQ(ierr); ierr = PetscViewerASCIIPrintf(viewer," Error estimation: %s\n",TSGLLEErrorDirections[gl->error_direction]);CHKERRQ(ierr); ierr = PetscViewerASCIIPrintf(viewer," Extrapolation: %s\n",gl->extrapolate ? "yes" : "no");CHKERRQ(ierr); ierr = PetscViewerASCIIPrintf(viewer," Acceptance test: %s\n",gl->accept_name[0] ? gl->accept_name : "(not yet set)");CHKERRQ(ierr); ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); ierr = TSGLLEAdaptView(gl->adapt,viewer);CHKERRQ(ierr); ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); ierr = PetscViewerASCIIPrintf(viewer," type: %s\n",gl->type_name[0] ? gl->type_name : "(not yet set)");CHKERRQ(ierr); ierr = PetscViewerASCIIPrintf(viewer,"Schemes within family (%d):\n",gl->nschemes);CHKERRQ(ierr); details = PETSC_FALSE; ierr = PetscOptionsGetBool(((PetscObject)ts)->options,((PetscObject)ts)->prefix,"-ts_gl_view_detailed",&details,NULL);CHKERRQ(ierr); ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); for (i=0; inschemes; i++) { ierr = TSGLLESchemeView(gl->schemes[i],details,viewer);CHKERRQ(ierr); } if (gl->View) { ierr = (*gl->View)(gl,viewer);CHKERRQ(ierr); } ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); } ierr = SNESView(ts->snes,viewer);CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "TSGLLERegister" /*@C TSGLLERegister - adds a TSGLLE implementation Not Collective Input Parameters: + name_scheme - name of user-defined general linear scheme - routine_create - routine to create method context Notes: TSGLLERegister() may be called multiple times to add several user-defined families. Sample usage: .vb TSGLLERegister("my_scheme",MySchemeCreate); .ve Then, your scheme can be chosen with the procedural interface via $ TSGLLESetType(ts,"my_scheme") or at runtime via the option $ -ts_gl_type my_scheme Level: advanced .keywords: TSGLLE, register .seealso: TSGLLERegisterAll() @*/ PetscErrorCode TSGLLERegister(const char sname[],PetscErrorCode (*function)(TS)) { PetscErrorCode ierr; PetscFunctionBegin; ierr = PetscFunctionListAdd(&TSGLLEList,sname,function);CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "TSGLLEAcceptRegister" /*@C TSGLLEAcceptRegister - adds a TSGLLE acceptance scheme Not Collective Input Parameters: + name_scheme - name of user-defined acceptance scheme - routine_create - routine to create method context Notes: TSGLLEAcceptRegister() may be called multiple times to add several user-defined families. Sample usage: .vb TSGLLEAcceptRegister("my_scheme",MySchemeCreate); .ve Then, your scheme can be chosen with the procedural interface via $ TSGLLESetAcceptType(ts,"my_scheme") or at runtime via the option $ -ts_gl_accept_type my_scheme Level: advanced .keywords: TSGLLE, TSGLLEAcceptType, register .seealso: TSGLLERegisterAll() @*/ PetscErrorCode TSGLLEAcceptRegister(const char sname[],TSGLLEAcceptFunction function) { PetscErrorCode ierr; PetscFunctionBegin; ierr = PetscFunctionListAdd(&TSGLLEAcceptList,sname,function);CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "TSGLLERegisterAll" /*@C TSGLLERegisterAll - Registers all of the general linear methods in TSGLLE Not Collective Level: advanced .keywords: TS, TSGLLE, register, all .seealso: TSGLLERegisterDestroy() @*/ PetscErrorCode TSGLLERegisterAll(void) { PetscErrorCode ierr; PetscFunctionBegin; if (TSGLLERegisterAllCalled) PetscFunctionReturn(0); TSGLLERegisterAllCalled = PETSC_TRUE; ierr = TSGLLERegister(TSGLLE_IRKS, TSGLLECreate_IRKS);CHKERRQ(ierr); ierr = TSGLLEAcceptRegister(TSGLLEACCEPT_ALWAYS,TSGLLEAccept_Always);CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "TSGLLEInitializePackage" /*@C TSGLLEInitializePackage - This function initializes everything in the TSGLLE package. It is called from PetscDLLibraryRegister() when using dynamic libraries, and on the first call to TSCreate_GLLE() when using static libraries. Level: developer .keywords: TS, TSGLLE, initialize, package .seealso: PetscInitialize() @*/ PetscErrorCode TSGLLEInitializePackage(void) { PetscErrorCode ierr; PetscFunctionBegin; if (TSGLLEPackageInitialized) PetscFunctionReturn(0); TSGLLEPackageInitialized = PETSC_TRUE; ierr = TSGLLERegisterAll();CHKERRQ(ierr); ierr = PetscRegisterFinalize(TSGLLEFinalizePackage);CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "TSGLLEFinalizePackage" /*@C TSGLLEFinalizePackage - This function destroys everything in the TSGLLE package. It is called from PetscFinalize(). Level: developer .keywords: Petsc, destroy, package .seealso: PetscFinalize() @*/ PetscErrorCode TSGLLEFinalizePackage(void) { PetscErrorCode ierr; PetscFunctionBegin; ierr = PetscFunctionListDestroy(&TSGLLEList);CHKERRQ(ierr); ierr = PetscFunctionListDestroy(&TSGLLEAcceptList);CHKERRQ(ierr); TSGLLEPackageInitialized = PETSC_FALSE; TSGLLERegisterAllCalled = PETSC_FALSE; PetscFunctionReturn(0); } /* ------------------------------------------------------------ */ /*MC TSGLLE - DAE solver using implicit General Linear methods These methods contain Runge-Kutta and multistep schemes as special cases. These special cases have some fundamental limitations. For example, diagonally implicit Runge-Kutta cannot have stage order greater than 1 which limits their applicability to very stiff systems. Meanwhile, multistep methods cannot be A-stable for order greater than 2 and BDF are not 0-stable for order greater than 6. GL methods can be A- and L-stable with arbitrarily high stage order and reliable error estimates for both 1 and 2 orders higher to facilitate adaptive step sizes and adaptive order schemes. All this is possible while preserving a singly diagonally implicit structure. Options database keys: + -ts_gl_type - the class of general linear method (irks) . -ts_gl_rtol - relative error . -ts_gl_atol - absolute error . -ts_gl_min_order

- minimum order method to consider (default=1) . -ts_gl_max_order

- maximum order method to consider (default=3) . -ts_gl_start_order

- order of starting method (default=1) . -ts_gl_complete - method to use for completing the step (rescale-and-modify or rescale) - -ts_adapt_type - adaptive controller to use (none step both) Notes: This integrator can be applied to DAE. Diagonally implicit general linear (DIGL) methods are a generalization of diagonally implicit Runge-Kutta (DIRK). They are represented by the tableau .vb A | U ------- B | V .ve combined with a vector c of abscissa. "Diagonally implicit" means that A is lower triangular. A step of the general method reads .vb [ Y ] = [A U] [ Y' ] [X^k] = [B V] [X^{k-1}] .ve where Y is the multivector of stage values, Y' is the multivector of stage derivatives, X^k is the Nordsieck vector of the solution at step k. The Nordsieck vector consists of the first r moments of the solution, given by .vb X = [x_0,x_1,...,x_{r-1}] = [x, h x', h^2 x'', ..., h^{r-1} x^{(r-1)} ] .ve If A is lower triangular, we can solve the stages (Y,Y') sequentially .vb 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} .ve and then construct the pieces to carry to the next step .vb 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} .ve Note that when the equations are cast in implicit form, we are using the stage equation to define y'_i in terms of y_i and known stuff (y_j for jdata = (void*)gl; ts->ops->reset = TSReset_GLLE; ts->ops->destroy = TSDestroy_GLLE; ts->ops->view = TSView_GLLE; ts->ops->setup = TSSetUp_GLLE; ts->ops->solve = TSSolve_GLLE; ts->ops->setfromoptions = TSSetFromOptions_GLLE; ts->ops->snesfunction = SNESTSFormFunction_GLLE; ts->ops->snesjacobian = SNESTSFormJacobian_GLLE; gl->max_step_rejections = 1; gl->min_order = 1; gl->max_order = 3; gl->start_order = 1; gl->current_scheme = -1; gl->extrapolate = PETSC_FALSE; gl->wrms_atol = 1e-8; gl->wrms_rtol = 1e-5; ierr = PetscObjectComposeFunction((PetscObject)ts,"TSGLLESetType_C", &TSGLLESetType_GLLE);CHKERRQ(ierr); ierr = PetscObjectComposeFunction((PetscObject)ts,"TSGLLESetAcceptType_C",&TSGLLESetAcceptType_GLLE);CHKERRQ(ierr); ierr = PetscObjectComposeFunction((PetscObject)ts,"TSGLLEGetAdapt_C", &TSGLLEGetAdapt_GLLE);CHKERRQ(ierr); PetscFunctionReturn(0); }