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