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