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 ts->total_steps++; 1001 1002 ierr = TSPostEvaluate(ts);CHKERRQ(ierr); 1003 ierr = TSPostStep(ts);CHKERRQ(ierr); 1004 ierr = TSMonitor(ts,ts->steps,ts->ptime,ts->vec_sol);CHKERRQ(ierr); 1005 1006 gl->current_scheme = next_scheme; 1007 ts->time_step = next_h; 1008 } 1009 PetscFunctionReturn(0); 1010 } 1011 1012 /*------------------------------------------------------------*/ 1013 1014 static PetscErrorCode TSReset_GLLE(TS ts) 1015 { 1016 TS_GLLE *gl = (TS_GLLE*)ts->data; 1017 PetscInt max_r,max_s; 1018 PetscErrorCode ierr; 1019 1020 PetscFunctionBegin; 1021 if (gl->setupcalled) { 1022 ierr = TSGLLEGetMaxSizes(ts,&max_r,&max_s);CHKERRQ(ierr); 1023 ierr = VecDestroyVecs(max_r,&gl->Xold);CHKERRQ(ierr); 1024 ierr = VecDestroyVecs(max_r,&gl->X);CHKERRQ(ierr); 1025 ierr = VecDestroyVecs(max_s,&gl->Ydot);CHKERRQ(ierr); 1026 ierr = VecDestroyVecs(3,&gl->himom);CHKERRQ(ierr); 1027 ierr = VecDestroy(&gl->W);CHKERRQ(ierr); 1028 ierr = VecDestroy(&gl->Y);CHKERRQ(ierr); 1029 ierr = VecDestroy(&gl->Z);CHKERRQ(ierr); 1030 } 1031 gl->setupcalled = PETSC_FALSE; 1032 PetscFunctionReturn(0); 1033 } 1034 1035 static PetscErrorCode TSDestroy_GLLE(TS ts) 1036 { 1037 TS_GLLE *gl = (TS_GLLE*)ts->data; 1038 PetscErrorCode ierr; 1039 1040 PetscFunctionBegin; 1041 ierr = TSReset_GLLE(ts);CHKERRQ(ierr); 1042 if (gl->adapt) {ierr = TSGLLEAdaptDestroy(&gl->adapt);CHKERRQ(ierr);} 1043 if (gl->Destroy) {ierr = (*gl->Destroy)(gl);CHKERRQ(ierr);} 1044 ierr = PetscFree(ts->data);CHKERRQ(ierr); 1045 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSGLLESetType_C",NULL);CHKERRQ(ierr); 1046 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSGLLESetAcceptType_C",NULL);CHKERRQ(ierr); 1047 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSGLLEGetAdapt_C",NULL);CHKERRQ(ierr); 1048 PetscFunctionReturn(0); 1049 } 1050 1051 /* 1052 This defines the nonlinear equation that is to be solved with SNES 1053 g(x) = f(t,x,z+shift*x) = 0 1054 */ 1055 static PetscErrorCode SNESTSFormFunction_GLLE(SNES snes,Vec x,Vec f,TS ts) 1056 { 1057 TS_GLLE *gl = (TS_GLLE*)ts->data; 1058 PetscErrorCode ierr; 1059 Vec Z,Ydot; 1060 DM dm,dmsave; 1061 1062 PetscFunctionBegin; 1063 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 1064 ierr = TSGLLEGetVecs(ts,dm,&Z,&Ydot);CHKERRQ(ierr); 1065 ierr = VecWAXPY(Ydot,gl->scoeff/ts->time_step,x,Z);CHKERRQ(ierr); 1066 dmsave = ts->dm; 1067 ts->dm = dm; 1068 ierr = TSComputeIFunction(ts,gl->stage_time,x,Ydot,f,PETSC_FALSE);CHKERRQ(ierr); 1069 ts->dm = dmsave; 1070 ierr = TSGLLERestoreVecs(ts,dm,&Z,&Ydot);CHKERRQ(ierr); 1071 PetscFunctionReturn(0); 1072 } 1073 1074 static PetscErrorCode SNESTSFormJacobian_GLLE(SNES snes,Vec x,Mat A,Mat B,TS ts) 1075 { 1076 TS_GLLE *gl = (TS_GLLE*)ts->data; 1077 PetscErrorCode ierr; 1078 Vec Z,Ydot; 1079 DM dm,dmsave; 1080 1081 PetscFunctionBegin; 1082 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 1083 ierr = TSGLLEGetVecs(ts,dm,&Z,&Ydot);CHKERRQ(ierr); 1084 dmsave = ts->dm; 1085 ts->dm = dm; 1086 /* gl->Xdot will have already been computed in SNESTSFormFunction_GLLE */ 1087 ierr = TSComputeIJacobian(ts,gl->stage_time,x,gl->Ydot[gl->stage],gl->scoeff/ts->time_step,A,B,PETSC_FALSE);CHKERRQ(ierr); 1088 ts->dm = dmsave; 1089 ierr = TSGLLERestoreVecs(ts,dm,&Z,&Ydot);CHKERRQ(ierr); 1090 PetscFunctionReturn(0); 1091 } 1092 1093 1094 static PetscErrorCode TSSetUp_GLLE(TS ts) 1095 { 1096 TS_GLLE *gl = (TS_GLLE*)ts->data; 1097 PetscInt max_r,max_s; 1098 PetscErrorCode ierr; 1099 DM dm; 1100 1101 PetscFunctionBegin; 1102 gl->setupcalled = PETSC_TRUE; 1103 ierr = TSGLLEGetMaxSizes(ts,&max_r,&max_s);CHKERRQ(ierr); 1104 ierr = VecDuplicateVecs(ts->vec_sol,max_r,&gl->X);CHKERRQ(ierr); 1105 ierr = VecDuplicateVecs(ts->vec_sol,max_r,&gl->Xold);CHKERRQ(ierr); 1106 ierr = VecDuplicateVecs(ts->vec_sol,max_s,&gl->Ydot);CHKERRQ(ierr); 1107 ierr = VecDuplicateVecs(ts->vec_sol,3,&gl->himom);CHKERRQ(ierr); 1108 ierr = VecDuplicate(ts->vec_sol,&gl->W);CHKERRQ(ierr); 1109 ierr = VecDuplicate(ts->vec_sol,&gl->Y);CHKERRQ(ierr); 1110 ierr = VecDuplicate(ts->vec_sol,&gl->Z);CHKERRQ(ierr); 1111 1112 /* Default acceptance tests and adaptivity */ 1113 if (!gl->Accept) {ierr = TSGLLESetAcceptType(ts,TSGLLEACCEPT_ALWAYS);CHKERRQ(ierr);} 1114 if (!gl->adapt) {ierr = TSGLLEGetAdapt(ts,&gl->adapt);CHKERRQ(ierr);} 1115 1116 if (gl->current_scheme < 0) { 1117 PetscInt i; 1118 for (i=0;; i++) { 1119 if (gl->schemes[i]->p == gl->start_order) break; 1120 if (i+1 == gl->nschemes) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"No schemes available with requested start order %d",i); 1121 } 1122 gl->current_scheme = i; 1123 } 1124 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 1125 if (dm) { 1126 ierr = DMCoarsenHookAdd(dm,DMCoarsenHook_TSGLLE,DMRestrictHook_TSGLLE,ts);CHKERRQ(ierr); 1127 ierr = DMSubDomainHookAdd(dm,DMSubDomainHook_TSGLLE,DMSubDomainRestrictHook_TSGLLE,ts);CHKERRQ(ierr); 1128 } 1129 PetscFunctionReturn(0); 1130 } 1131 /*------------------------------------------------------------*/ 1132 1133 static PetscErrorCode TSSetFromOptions_GLLE(PetscOptionItems *PetscOptionsObject,TS ts) 1134 { 1135 TS_GLLE *gl = (TS_GLLE*)ts->data; 1136 char tname[256] = TSGLLE_IRKS,completef[256] = "rescale-and-modify"; 1137 PetscErrorCode ierr; 1138 1139 PetscFunctionBegin; 1140 ierr = PetscOptionsHead(PetscOptionsObject,"General Linear ODE solver options");CHKERRQ(ierr); 1141 { 1142 PetscBool flg; 1143 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); 1144 if (flg || !gl->type_name[0]) { 1145 ierr = TSGLLESetType(ts,tname);CHKERRQ(ierr); 1146 } 1147 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); 1148 ierr = PetscOptionsInt("-ts_gl_max_order","Maximum order to try","TSGLLESetMaxOrder",gl->max_order,&gl->max_order,NULL);CHKERRQ(ierr); 1149 ierr = PetscOptionsInt("-ts_gl_min_order","Minimum order to try","TSGLLESetMinOrder",gl->min_order,&gl->min_order,NULL);CHKERRQ(ierr); 1150 ierr = PetscOptionsInt("-ts_gl_start_order","Initial order to try","TSGLLESetMinOrder",gl->start_order,&gl->start_order,NULL);CHKERRQ(ierr); 1151 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); 1152 ierr = PetscOptionsBool("-ts_gl_extrapolate","Extrapolate stage solution from previous solution (sometimes unstable)","TSGLLESetExtrapolate",gl->extrapolate,&gl->extrapolate,NULL);CHKERRQ(ierr); 1153 ierr = PetscOptionsReal("-ts_gl_atol","Absolute tolerance","TSGLLESetTolerances",gl->wrms_atol,&gl->wrms_atol,NULL);CHKERRQ(ierr); 1154 ierr = PetscOptionsReal("-ts_gl_rtol","Relative tolerance","TSGLLESetTolerances",gl->wrms_rtol,&gl->wrms_rtol,NULL);CHKERRQ(ierr); 1155 ierr = PetscOptionsString("-ts_gl_complete","Method to use for completing the step","none",completef,completef,sizeof(completef),&flg);CHKERRQ(ierr); 1156 if (flg) { 1157 PetscBool match1,match2; 1158 ierr = PetscStrcmp(completef,"rescale",&match1);CHKERRQ(ierr); 1159 ierr = PetscStrcmp(completef,"rescale-and-modify",&match2);CHKERRQ(ierr); 1160 if (match1) gl->CompleteStep = TSGLLECompleteStep_Rescale; 1161 else if (match2) gl->CompleteStep = TSGLLECompleteStep_RescaleAndModify; 1162 else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"%s",completef); 1163 } 1164 { 1165 char type[256] = TSGLLEACCEPT_ALWAYS; 1166 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); 1167 if (flg || !gl->accept_name[0]) { 1168 ierr = TSGLLESetAcceptType(ts,type);CHKERRQ(ierr); 1169 } 1170 } 1171 { 1172 TSGLLEAdapt adapt; 1173 ierr = TSGLLEGetAdapt(ts,&adapt);CHKERRQ(ierr); 1174 ierr = TSGLLEAdaptSetFromOptions(PetscOptionsObject,adapt);CHKERRQ(ierr); 1175 } 1176 } 1177 ierr = PetscOptionsTail();CHKERRQ(ierr); 1178 PetscFunctionReturn(0); 1179 } 1180 1181 static PetscErrorCode TSView_GLLE(TS ts,PetscViewer viewer) 1182 { 1183 TS_GLLE *gl = (TS_GLLE*)ts->data; 1184 PetscInt i; 1185 PetscBool iascii,details; 1186 PetscErrorCode ierr; 1187 1188 PetscFunctionBegin; 1189 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 1190 if (iascii) { 1191 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); 1192 ierr = PetscViewerASCIIPrintf(viewer," Error estimation: %s\n",TSGLLEErrorDirections[gl->error_direction]);CHKERRQ(ierr); 1193 ierr = PetscViewerASCIIPrintf(viewer," Extrapolation: %s\n",gl->extrapolate ? "yes" : "no");CHKERRQ(ierr); 1194 ierr = PetscViewerASCIIPrintf(viewer," Acceptance test: %s\n",gl->accept_name[0] ? gl->accept_name : "(not yet set)");CHKERRQ(ierr); 1195 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 1196 ierr = TSGLLEAdaptView(gl->adapt,viewer);CHKERRQ(ierr); 1197 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 1198 ierr = PetscViewerASCIIPrintf(viewer," type: %s\n",gl->type_name[0] ? gl->type_name : "(not yet set)");CHKERRQ(ierr); 1199 ierr = PetscViewerASCIIPrintf(viewer,"Schemes within family (%d):\n",gl->nschemes);CHKERRQ(ierr); 1200 details = PETSC_FALSE; 1201 ierr = PetscOptionsGetBool(((PetscObject)ts)->options,((PetscObject)ts)->prefix,"-ts_gl_view_detailed",&details,NULL);CHKERRQ(ierr); 1202 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 1203 for (i=0; i<gl->nschemes; i++) { 1204 ierr = TSGLLESchemeView(gl->schemes[i],details,viewer);CHKERRQ(ierr); 1205 } 1206 if (gl->View) { 1207 ierr = (*gl->View)(gl,viewer);CHKERRQ(ierr); 1208 } 1209 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 1210 } 1211 PetscFunctionReturn(0); 1212 } 1213 1214 /*@C 1215 TSGLLERegister - adds a TSGLLE implementation 1216 1217 Not Collective 1218 1219 Input Parameters: 1220 + name_scheme - name of user-defined general linear scheme 1221 - routine_create - routine to create method context 1222 1223 Notes: 1224 TSGLLERegister() may be called multiple times to add several user-defined families. 1225 1226 Sample usage: 1227 .vb 1228 TSGLLERegister("my_scheme",MySchemeCreate); 1229 .ve 1230 1231 Then, your scheme can be chosen with the procedural interface via 1232 $ TSGLLESetType(ts,"my_scheme") 1233 or at runtime via the option 1234 $ -ts_gl_type my_scheme 1235 1236 Level: advanced 1237 1238 .keywords: TSGLLE, register 1239 1240 .seealso: TSGLLERegisterAll() 1241 @*/ 1242 PetscErrorCode TSGLLERegister(const char sname[],PetscErrorCode (*function)(TS)) 1243 { 1244 PetscErrorCode ierr; 1245 1246 PetscFunctionBegin; 1247 ierr = PetscFunctionListAdd(&TSGLLEList,sname,function);CHKERRQ(ierr); 1248 PetscFunctionReturn(0); 1249 } 1250 1251 /*@C 1252 TSGLLEAcceptRegister - adds a TSGLLE acceptance scheme 1253 1254 Not Collective 1255 1256 Input Parameters: 1257 + name_scheme - name of user-defined acceptance scheme 1258 - routine_create - routine to create method context 1259 1260 Notes: 1261 TSGLLEAcceptRegister() may be called multiple times to add several user-defined families. 1262 1263 Sample usage: 1264 .vb 1265 TSGLLEAcceptRegister("my_scheme",MySchemeCreate); 1266 .ve 1267 1268 Then, your scheme can be chosen with the procedural interface via 1269 $ TSGLLESetAcceptType(ts,"my_scheme") 1270 or at runtime via the option 1271 $ -ts_gl_accept_type my_scheme 1272 1273 Level: advanced 1274 1275 .keywords: TSGLLE, TSGLLEAcceptType, register 1276 1277 .seealso: TSGLLERegisterAll() 1278 @*/ 1279 PetscErrorCode TSGLLEAcceptRegister(const char sname[],TSGLLEAcceptFunction function) 1280 { 1281 PetscErrorCode ierr; 1282 1283 PetscFunctionBegin; 1284 ierr = PetscFunctionListAdd(&TSGLLEAcceptList,sname,function);CHKERRQ(ierr); 1285 PetscFunctionReturn(0); 1286 } 1287 1288 /*@C 1289 TSGLLERegisterAll - Registers all of the general linear methods in TSGLLE 1290 1291 Not Collective 1292 1293 Level: advanced 1294 1295 .keywords: TS, TSGLLE, register, all 1296 1297 .seealso: TSGLLERegisterDestroy() 1298 @*/ 1299 PetscErrorCode TSGLLERegisterAll(void) 1300 { 1301 PetscErrorCode ierr; 1302 1303 PetscFunctionBegin; 1304 if (TSGLLERegisterAllCalled) PetscFunctionReturn(0); 1305 TSGLLERegisterAllCalled = PETSC_TRUE; 1306 1307 ierr = TSGLLERegister(TSGLLE_IRKS, TSGLLECreate_IRKS);CHKERRQ(ierr); 1308 ierr = TSGLLEAcceptRegister(TSGLLEACCEPT_ALWAYS,TSGLLEAccept_Always);CHKERRQ(ierr); 1309 PetscFunctionReturn(0); 1310 } 1311 1312 /*@C 1313 TSGLLEInitializePackage - This function initializes everything in the TSGLLE package. It is called 1314 from PetscDLLibraryRegister() when using dynamic libraries, and on the first call to TSCreate_GLLE() 1315 when using static libraries. 1316 1317 Level: developer 1318 1319 .keywords: TS, TSGLLE, initialize, package 1320 .seealso: PetscInitialize() 1321 @*/ 1322 PetscErrorCode TSGLLEInitializePackage(void) 1323 { 1324 PetscErrorCode ierr; 1325 1326 PetscFunctionBegin; 1327 if (TSGLLEPackageInitialized) PetscFunctionReturn(0); 1328 TSGLLEPackageInitialized = PETSC_TRUE; 1329 ierr = TSGLLERegisterAll();CHKERRQ(ierr); 1330 ierr = PetscRegisterFinalize(TSGLLEFinalizePackage);CHKERRQ(ierr); 1331 PetscFunctionReturn(0); 1332 } 1333 1334 /*@C 1335 TSGLLEFinalizePackage - This function destroys everything in the TSGLLE package. It is 1336 called from PetscFinalize(). 1337 1338 Level: developer 1339 1340 .keywords: Petsc, destroy, package 1341 .seealso: PetscFinalize() 1342 @*/ 1343 PetscErrorCode TSGLLEFinalizePackage(void) 1344 { 1345 PetscErrorCode ierr; 1346 1347 PetscFunctionBegin; 1348 ierr = PetscFunctionListDestroy(&TSGLLEList);CHKERRQ(ierr); 1349 ierr = PetscFunctionListDestroy(&TSGLLEAcceptList);CHKERRQ(ierr); 1350 TSGLLEPackageInitialized = PETSC_FALSE; 1351 TSGLLERegisterAllCalled = PETSC_FALSE; 1352 PetscFunctionReturn(0); 1353 } 1354 1355 /* ------------------------------------------------------------ */ 1356 /*MC 1357 TSGLLE - DAE solver using implicit General Linear methods 1358 1359 These methods contain Runge-Kutta and multistep schemes as special cases. These special cases have some fundamental 1360 limitations. For example, diagonally implicit Runge-Kutta cannot have stage order greater than 1 which limits their 1361 applicability to very stiff systems. Meanwhile, multistep methods cannot be A-stable for order greater than 2 and BDF 1362 are not 0-stable for order greater than 6. GL methods can be A- and L-stable with arbitrarily high stage order and 1363 reliable error estimates for both 1 and 2 orders higher to facilitate adaptive step sizes and adaptive order schemes. 1364 All this is possible while preserving a singly diagonally implicit structure. 1365 1366 Options database keys: 1367 + -ts_gl_type <type> - the class of general linear method (irks) 1368 . -ts_gl_rtol <tol> - relative error 1369 . -ts_gl_atol <tol> - absolute error 1370 . -ts_gl_min_order <p> - minimum order method to consider (default=1) 1371 . -ts_gl_max_order <p> - maximum order method to consider (default=3) 1372 . -ts_gl_start_order <p> - order of starting method (default=1) 1373 . -ts_gl_complete <method> - method to use for completing the step (rescale-and-modify or rescale) 1374 - -ts_adapt_type <method> - adaptive controller to use (none step both) 1375 1376 Notes: 1377 This integrator can be applied to DAE. 1378 1379 Diagonally implicit general linear (DIGL) methods are a generalization of diagonally implicit Runge-Kutta (DIRK). 1380 They are represented by the tableau 1381 1382 .vb 1383 A | U 1384 ------- 1385 B | V 1386 .ve 1387 1388 combined with a vector c of abscissa. "Diagonally implicit" means that A is lower triangular. 1389 A step of the general method reads 1390 1391 .vb 1392 [ Y ] = [A U] [ Y' ] 1393 [X^k] = [B V] [X^{k-1}] 1394 .ve 1395 1396 where Y is the multivector of stage values, Y' is the multivector of stage derivatives, X^k is the Nordsieck vector of 1397 the solution at step k. The Nordsieck vector consists of the first r moments of the solution, given by 1398 1399 .vb 1400 X = [x_0,x_1,...,x_{r-1}] = [x, h x', h^2 x'', ..., h^{r-1} x^{(r-1)} ] 1401 .ve 1402 1403 If A is lower triangular, we can solve the stages (Y,Y') sequentially 1404 1405 .vb 1406 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} 1407 .ve 1408 1409 and then construct the pieces to carry to the next step 1410 1411 .vb 1412 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} 1413 .ve 1414 1415 Note that when the equations are cast in implicit form, we are using the stage equation to define y'_i 1416 in terms of y_i and known stuff (y_j for j<i and x_j for all j). 1417 1418 1419 Error estimation 1420 1421 At present, the most attractive GL methods for stiff problems are singly diagonally implicit schemes which posses 1422 Inherent Runge-Kutta Stability (IRKS). These methods have r=s, the number of items passed between steps is equal to 1423 the number of stages. The order and stage-order are one less than the number of stages. We use the error estimates 1424 in the 2007 paper which provide the following estimates 1425 1426 .vb 1427 h^{p+1} X^{(p+1)} = phi_0^T Y' + [0 psi_0^T] Xold 1428 h^{p+2} X^{(p+2)} = phi_1^T Y' + [0 psi_1^T] Xold 1429 h^{p+2} (dx'/dx) X^{(p+1)} = phi_2^T Y' + [0 psi_2^T] Xold 1430 .ve 1431 1432 These estimates are accurate to O(h^{p+3}). 1433 1434 Changing the step size 1435 1436 We use the generalized "rescale and modify" scheme, see equation (4.5) of the 2007 paper. 1437 1438 Level: beginner 1439 1440 References: 1441 + 1. - John Butcher and Z. Jackieweicz and W. Wright, On error propagation in general linear methods for 1442 ordinary differential equations, Journal of Complexity, Vol 23, 2007. 1443 - 2. - John Butcher, Numerical methods for ordinary differential equations, second edition, Wiley, 2009. 1444 1445 .seealso: TSCreate(), TS, TSSetType() 1446 1447 M*/ 1448 PETSC_EXTERN PetscErrorCode TSCreate_GLLE(TS ts) 1449 { 1450 TS_GLLE *gl; 1451 PetscErrorCode ierr; 1452 1453 PetscFunctionBegin; 1454 ierr = TSGLLEInitializePackage();CHKERRQ(ierr); 1455 1456 ierr = PetscNewLog(ts,&gl);CHKERRQ(ierr); 1457 ts->data = (void*)gl; 1458 1459 ts->ops->reset = TSReset_GLLE; 1460 ts->ops->destroy = TSDestroy_GLLE; 1461 ts->ops->view = TSView_GLLE; 1462 ts->ops->setup = TSSetUp_GLLE; 1463 ts->ops->solve = TSSolve_GLLE; 1464 ts->ops->setfromoptions = TSSetFromOptions_GLLE; 1465 ts->ops->snesfunction = SNESTSFormFunction_GLLE; 1466 ts->ops->snesjacobian = SNESTSFormJacobian_GLLE; 1467 1468 ts->usessnes = PETSC_TRUE; 1469 1470 1471 gl->max_step_rejections = 1; 1472 gl->min_order = 1; 1473 gl->max_order = 3; 1474 gl->start_order = 1; 1475 gl->current_scheme = -1; 1476 gl->extrapolate = PETSC_FALSE; 1477 1478 gl->wrms_atol = 1e-8; 1479 gl->wrms_rtol = 1e-5; 1480 1481 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSGLLESetType_C", &TSGLLESetType_GLLE);CHKERRQ(ierr); 1482 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSGLLESetAcceptType_C",&TSGLLESetAcceptType_GLLE);CHKERRQ(ierr); 1483 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSGLLEGetAdapt_C", &TSGLLEGetAdapt_GLLE);CHKERRQ(ierr); 1484 PetscFunctionReturn(0); 1485 } 1486