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