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