1 2 #include <../src/ts/impls/implicit/glle/glle.h> /*I "petscts.h" I*/ 3 #include <petscdm.h> 4 #include <petscblaslapack.h> 5 6 static const char *TSGLLEErrorDirections[] = {"FORWARD","BACKWARD","TSGLLEErrorDirection","TSGLLEERROR_",0}; 7 static PetscFunctionList TSGLLEList; 8 static PetscFunctionList TSGLLEAcceptList; 9 static PetscBool TSGLLEPackageInitialized; 10 static PetscBool TSGLLERegisterAllCalled; 11 12 /* This function is pure */ 13 static PetscScalar Factorial(PetscInt n) 14 { 15 PetscInt i; 16 if (n < 12) { /* Can compute with 32-bit integers */ 17 PetscInt f = 1; 18 for (i=2; i<=n; i++) f *= i; 19 return (PetscScalar)f; 20 } else { 21 PetscScalar f = 1.; 22 for (i=2; i<=n; i++) f *= (PetscScalar)i; 23 return f; 24 } 25 } 26 27 /* This function is pure */ 28 static PetscScalar CPowF(PetscScalar c,PetscInt p) 29 { 30 return PetscPowRealInt(PetscRealPart(c),p)/Factorial(p); 31 } 32 33 static PetscErrorCode TSGLLEGetVecs(TS ts,DM dm,Vec *Z,Vec *Ydotstage) 34 { 35 TS_GLLE *gl = (TS_GLLE*)ts->data; 36 PetscErrorCode ierr; 37 38 PetscFunctionBegin; 39 if (Z) { 40 if (dm && dm != ts->dm) { 41 ierr = DMGetNamedGlobalVector(dm,"TSGLLE_Z",Z);CHKERRQ(ierr); 42 } else *Z = gl->Z; 43 } 44 if (Ydotstage) { 45 if (dm && dm != ts->dm) { 46 ierr = DMGetNamedGlobalVector(dm,"TSGLLE_Ydot",Ydotstage);CHKERRQ(ierr); 47 } else *Ydotstage = gl->Ydot[gl->stage]; 48 } 49 PetscFunctionReturn(0); 50 } 51 52 53 static PetscErrorCode TSGLLERestoreVecs(TS ts,DM dm,Vec *Z,Vec *Ydotstage) 54 { 55 PetscErrorCode ierr; 56 57 PetscFunctionBegin; 58 if (Z) { 59 if (dm && dm != ts->dm) { 60 ierr = DMRestoreNamedGlobalVector(dm,"TSGLLE_Z",Z);CHKERRQ(ierr); 61 } 62 } 63 if (Ydotstage) { 64 65 if (dm && dm != ts->dm) { 66 ierr = DMRestoreNamedGlobalVector(dm,"TSGLLE_Ydot",Ydotstage);CHKERRQ(ierr); 67 } 68 } 69 PetscFunctionReturn(0); 70 } 71 72 static PetscErrorCode DMCoarsenHook_TSGLLE(DM fine,DM coarse,void *ctx) 73 { 74 PetscFunctionBegin; 75 PetscFunctionReturn(0); 76 } 77 78 static PetscErrorCode DMRestrictHook_TSGLLE(DM fine,Mat restrct,Vec rscale,Mat inject,DM coarse,void *ctx) 79 { 80 TS ts = (TS)ctx; 81 PetscErrorCode ierr; 82 Vec Ydot,Ydot_c; 83 84 PetscFunctionBegin; 85 ierr = TSGLLEGetVecs(ts,fine,NULL,&Ydot);CHKERRQ(ierr); 86 ierr = TSGLLEGetVecs(ts,coarse,NULL,&Ydot_c);CHKERRQ(ierr); 87 ierr = MatRestrict(restrct,Ydot,Ydot_c);CHKERRQ(ierr); 88 ierr = VecPointwiseMult(Ydot_c,rscale,Ydot_c);CHKERRQ(ierr); 89 ierr = TSGLLERestoreVecs(ts,fine,NULL,&Ydot);CHKERRQ(ierr); 90 ierr = TSGLLERestoreVecs(ts,coarse,NULL,&Ydot_c);CHKERRQ(ierr); 91 PetscFunctionReturn(0); 92 } 93 94 static PetscErrorCode DMSubDomainHook_TSGLLE(DM dm,DM subdm,void *ctx) 95 { 96 PetscFunctionBegin; 97 PetscFunctionReturn(0); 98 } 99 100 static PetscErrorCode DMSubDomainRestrictHook_TSGLLE(DM dm,VecScatter gscat, VecScatter lscat,DM subdm,void *ctx) 101 { 102 TS ts = (TS)ctx; 103 PetscErrorCode ierr; 104 Vec Ydot,Ydot_s; 105 106 PetscFunctionBegin; 107 ierr = TSGLLEGetVecs(ts,dm,NULL,&Ydot);CHKERRQ(ierr); 108 ierr = TSGLLEGetVecs(ts,subdm,NULL,&Ydot_s);CHKERRQ(ierr); 109 110 ierr = VecScatterBegin(gscat,Ydot,Ydot_s,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 111 ierr = VecScatterEnd(gscat,Ydot,Ydot_s,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 112 113 ierr = TSGLLERestoreVecs(ts,dm,NULL,&Ydot);CHKERRQ(ierr); 114 ierr = TSGLLERestoreVecs(ts,subdm,NULL,&Ydot_s);CHKERRQ(ierr); 115 PetscFunctionReturn(0); 116 } 117 118 static PetscErrorCode TSGLLESchemeCreate(PetscInt p,PetscInt q,PetscInt r,PetscInt s,const PetscScalar *c, 119 const PetscScalar *a,const PetscScalar *b,const PetscScalar *u,const PetscScalar *v,TSGLLEScheme *inscheme) 120 { 121 TSGLLEScheme scheme; 122 PetscInt j; 123 PetscErrorCode ierr; 124 125 PetscFunctionBegin; 126 if (p < 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Scheme order must be positive"); 127 if (r < 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"At least one item must be carried between steps"); 128 if (s < 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"At least one stage is required"); 129 PetscValidPointer(inscheme,4); 130 *inscheme = 0; 131 ierr = PetscNew(&scheme);CHKERRQ(ierr); 132 scheme->p = p; 133 scheme->q = q; 134 scheme->r = r; 135 scheme->s = s; 136 137 ierr = PetscMalloc5(s,&scheme->c,s*s,&scheme->a,r*s,&scheme->b,r*s,&scheme->u,r*r,&scheme->v);CHKERRQ(ierr); 138 ierr = PetscArraycpy(scheme->c,c,s);CHKERRQ(ierr); 139 for (j=0; j<s*s; j++) scheme->a[j] = (PetscAbsScalar(a[j]) < 1e-12) ? 0 : a[j]; 140 for (j=0; j<r*s; j++) scheme->b[j] = (PetscAbsScalar(b[j]) < 1e-12) ? 0 : b[j]; 141 for (j=0; j<s*r; j++) scheme->u[j] = (PetscAbsScalar(u[j]) < 1e-12) ? 0 : u[j]; 142 for (j=0; j<r*r; j++) scheme->v[j] = (PetscAbsScalar(v[j]) < 1e-12) ? 0 : v[j]; 143 144 ierr = PetscMalloc6(r,&scheme->alpha,r,&scheme->beta,r,&scheme->gamma,3*s,&scheme->phi,3*r,&scheme->psi,r,&scheme->stage_error);CHKERRQ(ierr); 145 { 146 PetscInt i,j,k,ss=s+2; 147 PetscBLASInt m,n,one=1,*ipiv,lwork=4*((s+3)*3+3),info,ldb; 148 PetscReal rcond,*sing,*workreal; 149 PetscScalar *ImV,*H,*bmat,*workscalar,*c=scheme->c,*a=scheme->a,*b=scheme->b,*u=scheme->u,*v=scheme->v; 150 #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 @*/ 624 PetscErrorCode TSGLLESetType(TS ts,TSGLLEType type) 625 { 626 PetscErrorCode ierr; 627 628 PetscFunctionBegin; 629 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 630 PetscValidCharPointer(type,2); 631 ierr = PetscTryMethod(ts,"TSGLLESetType_C",(TS,TSGLLEType),(ts,type));CHKERRQ(ierr); 632 PetscFunctionReturn(0); 633 } 634 635 /*@C 636 TSGLLESetAcceptType - sets the acceptance test 637 638 Time integrators that need to control error must have the option to reject a time step based on local error 639 estimates. This function allows different schemes to be set. 640 641 Logically Collective on TS 642 643 Input Parameters: 644 + ts - the TS context 645 - type - the type 646 647 Options Database Key: 648 . -ts_gl_accept_type <type> - sets the method used to determine whether to accept or reject a step 649 650 Level: intermediate 651 652 .seealso: TS, TSGLLE, TSGLLEAcceptRegister(), TSGLLEAdapt, set type 653 @*/ 654 PetscErrorCode TSGLLESetAcceptType(TS ts,TSGLLEAcceptType type) 655 { 656 PetscErrorCode ierr; 657 658 PetscFunctionBegin; 659 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 660 PetscValidCharPointer(type,2); 661 ierr = PetscTryMethod(ts,"TSGLLESetAcceptType_C",(TS,TSGLLEAcceptType),(ts,type));CHKERRQ(ierr); 662 PetscFunctionReturn(0); 663 } 664 665 /*@C 666 TSGLLEGetAdapt - gets the TSGLLEAdapt object from the TS 667 668 Not Collective 669 670 Input Parameter: 671 . ts - the TS context 672 673 Output Parameter: 674 . adapt - the TSGLLEAdapt context 675 676 Notes: 677 This allows the user set options on the TSGLLEAdapt object. Usually it is better to do this using the options 678 database, so this function is rarely needed. 679 680 Level: advanced 681 682 .seealso: TSGLLEAdapt, TSGLLEAdaptRegister() 683 @*/ 684 PetscErrorCode TSGLLEGetAdapt(TS ts,TSGLLEAdapt *adapt) 685 { 686 PetscErrorCode ierr; 687 688 PetscFunctionBegin; 689 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 690 PetscValidPointer(adapt,2); 691 ierr = PetscUseMethod(ts,"TSGLLEGetAdapt_C",(TS,TSGLLEAdapt*),(ts,adapt));CHKERRQ(ierr); 692 PetscFunctionReturn(0); 693 } 694 695 static PetscErrorCode TSGLLEAccept_Always(TS ts,PetscReal tleft,PetscReal h,const PetscReal enorms[],PetscBool *accept) 696 { 697 PetscFunctionBegin; 698 *accept = PETSC_TRUE; 699 PetscFunctionReturn(0); 700 } 701 702 static PetscErrorCode TSGLLEUpdateWRMS(TS ts) 703 { 704 TS_GLLE *gl = (TS_GLLE*)ts->data; 705 PetscErrorCode ierr; 706 PetscScalar *x,*w; 707 PetscInt n,i; 708 709 PetscFunctionBegin; 710 ierr = VecGetArray(gl->X[0],&x);CHKERRQ(ierr); 711 ierr = VecGetArray(gl->W,&w);CHKERRQ(ierr); 712 ierr = VecGetLocalSize(gl->W,&n);CHKERRQ(ierr); 713 for (i=0; i<n; i++) w[i] = 1./(gl->wrms_atol + gl->wrms_rtol*PetscAbsScalar(x[i])); 714 ierr = VecRestoreArray(gl->X[0],&x);CHKERRQ(ierr); 715 ierr = VecRestoreArray(gl->W,&w);CHKERRQ(ierr); 716 PetscFunctionReturn(0); 717 } 718 719 static PetscErrorCode TSGLLEVecNormWRMS(TS ts,Vec X,PetscReal *nrm) 720 { 721 TS_GLLE *gl = (TS_GLLE*)ts->data; 722 PetscErrorCode ierr; 723 PetscScalar *x,*w; 724 PetscReal sum = 0.0,gsum; 725 PetscInt n,N,i; 726 727 PetscFunctionBegin; 728 ierr = VecGetArray(X,&x);CHKERRQ(ierr); 729 ierr = VecGetArray(gl->W,&w);CHKERRQ(ierr); 730 ierr = VecGetLocalSize(gl->W,&n);CHKERRQ(ierr); 731 for (i=0; i<n; i++) sum += PetscAbsScalar(PetscSqr(x[i]*w[i])); 732 ierr = VecRestoreArray(X,&x);CHKERRQ(ierr); 733 ierr = VecRestoreArray(gl->W,&w);CHKERRQ(ierr); 734 ierr = MPIU_Allreduce(&sum,&gsum,1,MPIU_REAL,MPIU_SUM,PetscObjectComm((PetscObject)ts));CHKERRQ(ierr); 735 ierr = VecGetSize(gl->W,&N);CHKERRQ(ierr); 736 *nrm = PetscSqrtReal(gsum/(1.*N)); 737 PetscFunctionReturn(0); 738 } 739 740 static PetscErrorCode TSGLLESetType_GLLE(TS ts,TSGLLEType type) 741 { 742 PetscErrorCode ierr,(*r)(TS); 743 PetscBool same; 744 TS_GLLE *gl = (TS_GLLE*)ts->data; 745 746 PetscFunctionBegin; 747 if (gl->type_name[0]) { 748 ierr = PetscStrcmp(gl->type_name,type,&same);CHKERRQ(ierr); 749 if (same) PetscFunctionReturn(0); 750 ierr = (*gl->Destroy)(gl);CHKERRQ(ierr); 751 } 752 753 ierr = PetscFunctionListFind(TSGLLEList,type,&r);CHKERRQ(ierr); 754 if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown TSGLLE type \"%s\" given",type); 755 ierr = (*r)(ts);CHKERRQ(ierr); 756 ierr = PetscStrcpy(gl->type_name,type);CHKERRQ(ierr); 757 PetscFunctionReturn(0); 758 } 759 760 static PetscErrorCode TSGLLESetAcceptType_GLLE(TS ts,TSGLLEAcceptType type) 761 { 762 PetscErrorCode ierr; 763 TSGLLEAcceptFunction r; 764 TS_GLLE *gl = (TS_GLLE*)ts->data; 765 766 PetscFunctionBegin; 767 ierr = PetscFunctionListFind(TSGLLEAcceptList,type,&r);CHKERRQ(ierr); 768 if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown TSGLLEAccept type \"%s\" given",type); 769 gl->Accept = r; 770 ierr = PetscStrncpy(gl->accept_name,type,sizeof(gl->accept_name));CHKERRQ(ierr); 771 PetscFunctionReturn(0); 772 } 773 774 static PetscErrorCode TSGLLEGetAdapt_GLLE(TS ts,TSGLLEAdapt *adapt) 775 { 776 PetscErrorCode ierr; 777 TS_GLLE *gl = (TS_GLLE*)ts->data; 778 779 PetscFunctionBegin; 780 if (!gl->adapt) { 781 ierr = TSGLLEAdaptCreate(PetscObjectComm((PetscObject)ts),&gl->adapt);CHKERRQ(ierr); 782 ierr = PetscObjectIncrementTabLevel((PetscObject)gl->adapt,(PetscObject)ts,1);CHKERRQ(ierr); 783 ierr = PetscLogObjectParent((PetscObject)ts,(PetscObject)gl->adapt);CHKERRQ(ierr); 784 } 785 *adapt = gl->adapt; 786 PetscFunctionReturn(0); 787 } 788 789 static PetscErrorCode TSGLLEChooseNextScheme(TS ts,PetscReal h,const PetscReal hmnorm[],PetscInt *next_scheme,PetscReal *next_h,PetscBool *finish) 790 { 791 PetscErrorCode ierr; 792 TS_GLLE *gl = (TS_GLLE*)ts->data; 793 PetscInt i,n,cur_p,cur,next_sc,candidates[64],orders[64]; 794 PetscReal errors[64],costs[64],tleft; 795 796 PetscFunctionBegin; 797 cur = -1; 798 cur_p = gl->schemes[gl->current_scheme]->p; 799 tleft = ts->max_time - (ts->ptime + ts->time_step); 800 for (i=0,n=0; i<gl->nschemes; i++) { 801 TSGLLEScheme sc = gl->schemes[i]; 802 if (sc->p < gl->min_order || gl->max_order < sc->p) continue; 803 if (sc->p == cur_p - 1) errors[n] = PetscAbsScalar(sc->alpha[0])*hmnorm[0]; 804 else if (sc->p == cur_p) errors[n] = PetscAbsScalar(sc->alpha[0])*hmnorm[1]; 805 else if (sc->p == cur_p+1) errors[n] = PetscAbsScalar(sc->alpha[0])*(hmnorm[2]+hmnorm[3]); 806 else continue; 807 candidates[n] = i; 808 orders[n] = PetscMin(sc->p,sc->q); /* order of global truncation error */ 809 costs[n] = sc->s; /* estimate the cost as the number of stages */ 810 if (i == gl->current_scheme) cur = n; 811 n++; 812 } 813 if (cur < 0 || gl->nschemes <= cur) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Current scheme not found in scheme list"); 814 ierr = TSGLLEAdaptChoose(gl->adapt,n,orders,errors,costs,cur,h,tleft,&next_sc,next_h,finish);CHKERRQ(ierr); 815 *next_scheme = candidates[next_sc]; 816 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); 817 PetscFunctionReturn(0); 818 } 819 820 static PetscErrorCode TSGLLEGetMaxSizes(TS ts,PetscInt *max_r,PetscInt *max_s) 821 { 822 TS_GLLE *gl = (TS_GLLE*)ts->data; 823 824 PetscFunctionBegin; 825 *max_r = gl->schemes[gl->nschemes-1]->r; 826 *max_s = gl->schemes[gl->nschemes-1]->s; 827 PetscFunctionReturn(0); 828 } 829 830 static PetscErrorCode TSSolve_GLLE(TS ts) 831 { 832 TS_GLLE *gl = (TS_GLLE*)ts->data; 833 PetscInt i,k,its,lits,max_r,max_s; 834 PetscBool final_step,finish; 835 SNESConvergedReason snesreason; 836 PetscErrorCode ierr; 837 838 PetscFunctionBegin; 839 ierr = TSMonitor(ts,ts->steps,ts->ptime,ts->vec_sol);CHKERRQ(ierr); 840 841 ierr = TSGLLEGetMaxSizes(ts,&max_r,&max_s);CHKERRQ(ierr); 842 ierr = VecCopy(ts->vec_sol,gl->X[0]);CHKERRQ(ierr); 843 for (i=1; i<max_r; i++) { 844 ierr = VecZeroEntries(gl->X[i]);CHKERRQ(ierr); 845 } 846 ierr = TSGLLEUpdateWRMS(ts);CHKERRQ(ierr); 847 848 if (0) { 849 /* Find consistent initial data for DAE */ 850 gl->stage_time = ts->ptime + ts->time_step; 851 gl->scoeff = 1.; 852 gl->stage = 0; 853 854 ierr = VecCopy(ts->vec_sol,gl->Z);CHKERRQ(ierr); 855 ierr = VecCopy(ts->vec_sol,gl->Y);CHKERRQ(ierr); 856 ierr = SNESSolve(ts->snes,NULL,gl->Y);CHKERRQ(ierr); 857 ierr = SNESGetIterationNumber(ts->snes,&its);CHKERRQ(ierr); 858 ierr = SNESGetLinearSolveIterations(ts->snes,&lits);CHKERRQ(ierr); 859 ierr = SNESGetConvergedReason(ts->snes,&snesreason);CHKERRQ(ierr); 860 861 ts->snes_its += its; ts->ksp_its += lits; 862 if (snesreason < 0 && ts->max_snes_failures > 0 && ++ts->num_snes_failures >= ts->max_snes_failures) { 863 ts->reason = TS_DIVERGED_NONLINEAR_SOLVE; 864 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); 865 PetscFunctionReturn(0); 866 } 867 } 868 869 if (gl->current_scheme < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"A starting scheme has not been provided"); 870 871 for (k=0,final_step=PETSC_FALSE,finish=PETSC_FALSE; k<ts->max_steps && !finish; k++) { 872 PetscInt j,r,s,next_scheme = 0,rejections; 873 PetscReal h,hmnorm[4],enorm[3],next_h; 874 PetscBool accept; 875 const PetscScalar *c,*a,*u; 876 Vec *X,*Ydot,Y; 877 TSGLLEScheme scheme = gl->schemes[gl->current_scheme]; 878 879 r = scheme->r; s = scheme->s; 880 c = scheme->c; 881 a = scheme->a; u = scheme->u; 882 h = ts->time_step; 883 X = gl->X; Ydot = gl->Ydot; Y = gl->Y; 884 885 if (ts->ptime > ts->max_time) break; 886 887 /* 888 We only call PreStep at the start of each STEP, not each STAGE. This is because it is 889 possible to fail (have to restart a step) after multiple stages. 890 */ 891 ierr = TSPreStep(ts);CHKERRQ(ierr); 892 893 rejections = 0; 894 while (1) { 895 for (i=0; i<s; i++) { 896 PetscScalar shift; 897 gl->scoeff = 1./PetscRealPart(a[i*s+i]); 898 shift = gl->scoeff/ts->time_step; 899 gl->stage = i; 900 gl->stage_time = ts->ptime + PetscRealPart(c[i])*h; 901 902 /* 903 * Stage equation: Y = h A Y' + U X 904 * We assume that A is lower-triangular so that we can solve the stages (Y,Y') sequentially 905 * Build the affine vector z_i = -[1/(h a_ii)](h sum_j a_ij y'_j + sum_j u_ij x_j) 906 * Then y'_i = z + 1/(h a_ii) y_i 907 */ 908 ierr = VecZeroEntries(gl->Z);CHKERRQ(ierr); 909 for (j=0; j<r; j++) { 910 ierr = VecAXPY(gl->Z,-shift*u[i*r+j],X[j]);CHKERRQ(ierr); 911 } 912 for (j=0; j<i; j++) { 913 ierr = VecAXPY(gl->Z,-shift*h*a[i*s+j],Ydot[j]);CHKERRQ(ierr); 914 } 915 /* Note: Z is used within function evaluation, Ydot = Z + shift*Y */ 916 917 /* Compute an estimate of Y to start Newton iteration */ 918 if (gl->extrapolate) { 919 if (i==0) { 920 /* Linear extrapolation on the first stage */ 921 ierr = VecWAXPY(Y,c[i]*h,X[1],X[0]);CHKERRQ(ierr); 922 } else { 923 /* Linear extrapolation from the last stage */ 924 ierr = VecAXPY(Y,(c[i]-c[i-1])*h,Ydot[i-1]);CHKERRQ(ierr); 925 } 926 } else if (i==0) { /* Directly use solution from the last step, otherwise reuse the last stage (do nothing) */ 927 ierr = VecCopy(X[0],Y);CHKERRQ(ierr); 928 } 929 930 /* Solve this stage (Ydot[i] is computed during function evaluation) */ 931 ierr = SNESSolve(ts->snes,NULL,Y);CHKERRQ(ierr); 932 ierr = SNESGetIterationNumber(ts->snes,&its);CHKERRQ(ierr); 933 ierr = SNESGetLinearSolveIterations(ts->snes,&lits);CHKERRQ(ierr); 934 ierr = SNESGetConvergedReason(ts->snes,&snesreason);CHKERRQ(ierr); 935 ts->snes_its += its; ts->ksp_its += lits; 936 if (snesreason < 0 && ts->max_snes_failures > 0 && ++ts->num_snes_failures >= ts->max_snes_failures) { 937 ts->reason = TS_DIVERGED_NONLINEAR_SOLVE; 938 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); 939 PetscFunctionReturn(0); 940 } 941 } 942 943 gl->stage_time = ts->ptime + ts->time_step; 944 945 ierr = (*gl->EstimateHigherMoments)(scheme,h,Ydot,gl->X,gl->himom);CHKERRQ(ierr); 946 /* hmnorm[i] = h^{p+i}x^{(p+i)} with i=0,1,2; hmnorm[3] = h^{p+2}(dx'/dx) x^{(p+1)} */ 947 for (i=0; i<3; i++) { 948 ierr = TSGLLEVecNormWRMS(ts,gl->himom[i],&hmnorm[i+1]);CHKERRQ(ierr); 949 } 950 enorm[0] = PetscRealPart(scheme->alpha[0])*hmnorm[1]; 951 enorm[1] = PetscRealPart(scheme->beta[0]) *hmnorm[2]; 952 enorm[2] = PetscRealPart(scheme->gamma[0])*hmnorm[3]; 953 ierr = (*gl->Accept)(ts,ts->max_time-gl->stage_time,h,enorm,&accept);CHKERRQ(ierr); 954 if (accept) goto accepted; 955 rejections++; 956 ierr = PetscInfo3(ts,"Step %D (t=%g) not accepted, rejections=%D\n",k,gl->stage_time,rejections);CHKERRQ(ierr); 957 if (rejections > gl->max_step_rejections) break; 958 /* 959 There are lots of reasons why a step might be rejected, including solvers not converging and other factors that 960 TSGLLEChooseNextScheme does not support. Additionally, the error estimates may be very screwed up, so I'm not 961 convinced that it's safe to just compute a new error estimate using the same interface as the current adaptor 962 (the adaptor interface probably has to change). Here we make an arbitrary and naive choice. This assumes that 963 steps were written in Nordsieck form. The "correct" method would be to re-complete the previous time step with 964 the correct "next" step size. It is unclear to me whether the present ad-hoc method of rescaling X is stable. 965 */ 966 h *= 0.5; 967 for (i=1; i<scheme->r; i++) { 968 ierr = VecScale(X[i],PetscPowRealInt(0.5,i));CHKERRQ(ierr); 969 } 970 } 971 SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_CONV_FAILED,"Time step %D (t=%g) not accepted after %D failures\n",k,gl->stage_time,rejections); 972 973 accepted: 974 /* This term is not error, but it *would* be the leading term for a lower order method */ 975 ierr = TSGLLEVecNormWRMS(ts,gl->X[scheme->r-1],&hmnorm[0]);CHKERRQ(ierr); 976 /* Correct scaling so that these are equivalent to norms of the Nordsieck vectors */ 977 978 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); 979 if (!final_step) { 980 ierr = TSGLLEChooseNextScheme(ts,h,hmnorm,&next_scheme,&next_h,&final_step);CHKERRQ(ierr); 981 } else { 982 /* Dummy values to complete the current step in a consistent manner */ 983 next_scheme = gl->current_scheme; 984 next_h = h; 985 finish = PETSC_TRUE; 986 } 987 988 X = gl->Xold; 989 gl->Xold = gl->X; 990 gl->X = X; 991 ierr = (*gl->CompleteStep)(scheme,h,gl->schemes[next_scheme],next_h,Ydot,gl->Xold,gl->X);CHKERRQ(ierr); 992 993 ierr = TSGLLEUpdateWRMS(ts);CHKERRQ(ierr); 994 995 /* Post the solution for the user, we could avoid this copy with a small bit of cleverness */ 996 ierr = VecCopy(gl->X[0],ts->vec_sol);CHKERRQ(ierr); 997 ts->ptime += h; 998 ts->steps++; 999 1000 ierr = TSPostEvaluate(ts);CHKERRQ(ierr); 1001 ierr = TSPostStep(ts);CHKERRQ(ierr); 1002 ierr = TSMonitor(ts,ts->steps,ts->ptime,ts->vec_sol);CHKERRQ(ierr); 1003 1004 gl->current_scheme = next_scheme; 1005 ts->time_step = next_h; 1006 } 1007 PetscFunctionReturn(0); 1008 } 1009 1010 /*------------------------------------------------------------*/ 1011 1012 static PetscErrorCode TSReset_GLLE(TS ts) 1013 { 1014 TS_GLLE *gl = (TS_GLLE*)ts->data; 1015 PetscInt max_r,max_s; 1016 PetscErrorCode ierr; 1017 1018 PetscFunctionBegin; 1019 if (gl->setupcalled) { 1020 ierr = TSGLLEGetMaxSizes(ts,&max_r,&max_s);CHKERRQ(ierr); 1021 ierr = VecDestroyVecs(max_r,&gl->Xold);CHKERRQ(ierr); 1022 ierr = VecDestroyVecs(max_r,&gl->X);CHKERRQ(ierr); 1023 ierr = VecDestroyVecs(max_s,&gl->Ydot);CHKERRQ(ierr); 1024 ierr = VecDestroyVecs(3,&gl->himom);CHKERRQ(ierr); 1025 ierr = VecDestroy(&gl->W);CHKERRQ(ierr); 1026 ierr = VecDestroy(&gl->Y);CHKERRQ(ierr); 1027 ierr = VecDestroy(&gl->Z);CHKERRQ(ierr); 1028 } 1029 gl->setupcalled = PETSC_FALSE; 1030 PetscFunctionReturn(0); 1031 } 1032 1033 static PetscErrorCode TSDestroy_GLLE(TS ts) 1034 { 1035 TS_GLLE *gl = (TS_GLLE*)ts->data; 1036 PetscErrorCode ierr; 1037 1038 PetscFunctionBegin; 1039 ierr = TSReset_GLLE(ts);CHKERRQ(ierr); 1040 if (ts->dm) { 1041 ierr = DMCoarsenHookRemove(ts->dm,DMCoarsenHook_TSGLLE,DMRestrictHook_TSGLLE,ts);CHKERRQ(ierr); 1042 ierr = DMSubDomainHookRemove(ts->dm,DMSubDomainHook_TSGLLE,DMSubDomainRestrictHook_TSGLLE,ts);CHKERRQ(ierr); 1043 } 1044 if (gl->adapt) {ierr = TSGLLEAdaptDestroy(&gl->adapt);CHKERRQ(ierr);} 1045 if (gl->Destroy) {ierr = (*gl->Destroy)(gl);CHKERRQ(ierr);} 1046 ierr = PetscFree(ts->data);CHKERRQ(ierr); 1047 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSGLLESetType_C",NULL);CHKERRQ(ierr); 1048 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSGLLESetAcceptType_C",NULL);CHKERRQ(ierr); 1049 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSGLLEGetAdapt_C",NULL);CHKERRQ(ierr); 1050 PetscFunctionReturn(0); 1051 } 1052 1053 /* 1054 This defines the nonlinear equation that is to be solved with SNES 1055 g(x) = f(t,x,z+shift*x) = 0 1056 */ 1057 static PetscErrorCode SNESTSFormFunction_GLLE(SNES snes,Vec x,Vec f,TS ts) 1058 { 1059 TS_GLLE *gl = (TS_GLLE*)ts->data; 1060 PetscErrorCode ierr; 1061 Vec Z,Ydot; 1062 DM dm,dmsave; 1063 1064 PetscFunctionBegin; 1065 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 1066 ierr = TSGLLEGetVecs(ts,dm,&Z,&Ydot);CHKERRQ(ierr); 1067 ierr = VecWAXPY(Ydot,gl->scoeff/ts->time_step,x,Z);CHKERRQ(ierr); 1068 dmsave = ts->dm; 1069 ts->dm = dm; 1070 ierr = TSComputeIFunction(ts,gl->stage_time,x,Ydot,f,PETSC_FALSE);CHKERRQ(ierr); 1071 ts->dm = dmsave; 1072 ierr = TSGLLERestoreVecs(ts,dm,&Z,&Ydot);CHKERRQ(ierr); 1073 PetscFunctionReturn(0); 1074 } 1075 1076 static PetscErrorCode SNESTSFormJacobian_GLLE(SNES snes,Vec x,Mat A,Mat B,TS ts) 1077 { 1078 TS_GLLE *gl = (TS_GLLE*)ts->data; 1079 PetscErrorCode ierr; 1080 Vec Z,Ydot; 1081 DM dm,dmsave; 1082 1083 PetscFunctionBegin; 1084 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 1085 ierr = TSGLLEGetVecs(ts,dm,&Z,&Ydot);CHKERRQ(ierr); 1086 dmsave = ts->dm; 1087 ts->dm = dm; 1088 /* gl->Xdot will have already been computed in SNESTSFormFunction_GLLE */ 1089 ierr = TSComputeIJacobian(ts,gl->stage_time,x,gl->Ydot[gl->stage],gl->scoeff/ts->time_step,A,B,PETSC_FALSE);CHKERRQ(ierr); 1090 ts->dm = dmsave; 1091 ierr = TSGLLERestoreVecs(ts,dm,&Z,&Ydot);CHKERRQ(ierr); 1092 PetscFunctionReturn(0); 1093 } 1094 1095 1096 static PetscErrorCode TSSetUp_GLLE(TS ts) 1097 { 1098 TS_GLLE *gl = (TS_GLLE*)ts->data; 1099 PetscInt max_r,max_s; 1100 PetscErrorCode ierr; 1101 DM dm; 1102 1103 PetscFunctionBegin; 1104 gl->setupcalled = PETSC_TRUE; 1105 ierr = TSGLLEGetMaxSizes(ts,&max_r,&max_s);CHKERRQ(ierr); 1106 ierr = VecDuplicateVecs(ts->vec_sol,max_r,&gl->X);CHKERRQ(ierr); 1107 ierr = VecDuplicateVecs(ts->vec_sol,max_r,&gl->Xold);CHKERRQ(ierr); 1108 ierr = VecDuplicateVecs(ts->vec_sol,max_s,&gl->Ydot);CHKERRQ(ierr); 1109 ierr = VecDuplicateVecs(ts->vec_sol,3,&gl->himom);CHKERRQ(ierr); 1110 ierr = VecDuplicate(ts->vec_sol,&gl->W);CHKERRQ(ierr); 1111 ierr = VecDuplicate(ts->vec_sol,&gl->Y);CHKERRQ(ierr); 1112 ierr = VecDuplicate(ts->vec_sol,&gl->Z);CHKERRQ(ierr); 1113 1114 /* Default acceptance tests and adaptivity */ 1115 if (!gl->Accept) {ierr = TSGLLESetAcceptType(ts,TSGLLEACCEPT_ALWAYS);CHKERRQ(ierr);} 1116 if (!gl->adapt) {ierr = TSGLLEGetAdapt(ts,&gl->adapt);CHKERRQ(ierr);} 1117 1118 if (gl->current_scheme < 0) { 1119 PetscInt i; 1120 for (i=0;; i++) { 1121 if (gl->schemes[i]->p == gl->start_order) break; 1122 if (i+1 == gl->nschemes) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"No schemes available with requested start order %d",i); 1123 } 1124 gl->current_scheme = i; 1125 } 1126 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 1127 ierr = DMCoarsenHookAdd(dm,DMCoarsenHook_TSGLLE,DMRestrictHook_TSGLLE,ts);CHKERRQ(ierr); 1128 ierr = DMSubDomainHookAdd(dm,DMSubDomainHook_TSGLLE,DMSubDomainRestrictHook_TSGLLE,ts);CHKERRQ(ierr); 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 .seealso: TSGLLERegisterAll() 1239 @*/ 1240 PetscErrorCode TSGLLERegister(const char sname[],PetscErrorCode (*function)(TS)) 1241 { 1242 PetscErrorCode ierr; 1243 1244 PetscFunctionBegin; 1245 ierr = TSGLLEInitializePackage();CHKERRQ(ierr); 1246 ierr = PetscFunctionListAdd(&TSGLLEList,sname,function);CHKERRQ(ierr); 1247 PetscFunctionReturn(0); 1248 } 1249 1250 /*@C 1251 TSGLLEAcceptRegister - adds a TSGLLE acceptance scheme 1252 1253 Not Collective 1254 1255 Input Parameters: 1256 + name_scheme - name of user-defined acceptance scheme 1257 - routine_create - routine to create method context 1258 1259 Notes: 1260 TSGLLEAcceptRegister() may be called multiple times to add several user-defined families. 1261 1262 Sample usage: 1263 .vb 1264 TSGLLEAcceptRegister("my_scheme",MySchemeCreate); 1265 .ve 1266 1267 Then, your scheme can be chosen with the procedural interface via 1268 $ TSGLLESetAcceptType(ts,"my_scheme") 1269 or at runtime via the option 1270 $ -ts_gl_accept_type my_scheme 1271 1272 Level: advanced 1273 1274 .seealso: TSGLLERegisterAll() 1275 @*/ 1276 PetscErrorCode TSGLLEAcceptRegister(const char sname[],TSGLLEAcceptFunction function) 1277 { 1278 PetscErrorCode ierr; 1279 1280 PetscFunctionBegin; 1281 ierr = PetscFunctionListAdd(&TSGLLEAcceptList,sname,function);CHKERRQ(ierr); 1282 PetscFunctionReturn(0); 1283 } 1284 1285 /*@C 1286 TSGLLERegisterAll - Registers all of the general linear methods in TSGLLE 1287 1288 Not Collective 1289 1290 Level: advanced 1291 1292 .seealso: TSGLLERegisterDestroy() 1293 @*/ 1294 PetscErrorCode TSGLLERegisterAll(void) 1295 { 1296 PetscErrorCode ierr; 1297 1298 PetscFunctionBegin; 1299 if (TSGLLERegisterAllCalled) PetscFunctionReturn(0); 1300 TSGLLERegisterAllCalled = PETSC_TRUE; 1301 1302 ierr = TSGLLERegister(TSGLLE_IRKS, TSGLLECreate_IRKS);CHKERRQ(ierr); 1303 ierr = TSGLLEAcceptRegister(TSGLLEACCEPT_ALWAYS,TSGLLEAccept_Always);CHKERRQ(ierr); 1304 PetscFunctionReturn(0); 1305 } 1306 1307 /*@C 1308 TSGLLEInitializePackage - This function initializes everything in the TSGLLE package. It is called 1309 from TSInitializePackage(). 1310 1311 Level: developer 1312 1313 .seealso: PetscInitialize() 1314 @*/ 1315 PetscErrorCode TSGLLEInitializePackage(void) 1316 { 1317 PetscErrorCode ierr; 1318 1319 PetscFunctionBegin; 1320 if (TSGLLEPackageInitialized) PetscFunctionReturn(0); 1321 TSGLLEPackageInitialized = PETSC_TRUE; 1322 ierr = TSGLLERegisterAll();CHKERRQ(ierr); 1323 ierr = PetscRegisterFinalize(TSGLLEFinalizePackage);CHKERRQ(ierr); 1324 PetscFunctionReturn(0); 1325 } 1326 1327 /*@C 1328 TSGLLEFinalizePackage - This function destroys everything in the TSGLLE package. It is 1329 called from PetscFinalize(). 1330 1331 Level: developer 1332 1333 .seealso: PetscFinalize() 1334 @*/ 1335 PetscErrorCode TSGLLEFinalizePackage(void) 1336 { 1337 PetscErrorCode ierr; 1338 1339 PetscFunctionBegin; 1340 ierr = PetscFunctionListDestroy(&TSGLLEList);CHKERRQ(ierr); 1341 ierr = PetscFunctionListDestroy(&TSGLLEAcceptList);CHKERRQ(ierr); 1342 TSGLLEPackageInitialized = PETSC_FALSE; 1343 TSGLLERegisterAllCalled = PETSC_FALSE; 1344 PetscFunctionReturn(0); 1345 } 1346 1347 /* ------------------------------------------------------------ */ 1348 /*MC 1349 TSGLLE - DAE solver using implicit General Linear methods 1350 1351 These methods contain Runge-Kutta and multistep schemes as special cases. These special cases have some fundamental 1352 limitations. For example, diagonally implicit Runge-Kutta cannot have stage order greater than 1 which limits their 1353 applicability to very stiff systems. Meanwhile, multistep methods cannot be A-stable for order greater than 2 and BDF 1354 are not 0-stable for order greater than 6. GL methods can be A- and L-stable with arbitrarily high stage order and 1355 reliable error estimates for both 1 and 2 orders higher to facilitate adaptive step sizes and adaptive order schemes. 1356 All this is possible while preserving a singly diagonally implicit structure. 1357 1358 Options database keys: 1359 + -ts_gl_type <type> - the class of general linear method (irks) 1360 . -ts_gl_rtol <tol> - relative error 1361 . -ts_gl_atol <tol> - absolute error 1362 . -ts_gl_min_order <p> - minimum order method to consider (default=1) 1363 . -ts_gl_max_order <p> - maximum order method to consider (default=3) 1364 . -ts_gl_start_order <p> - order of starting method (default=1) 1365 . -ts_gl_complete <method> - method to use for completing the step (rescale-and-modify or rescale) 1366 - -ts_adapt_type <method> - adaptive controller to use (none step both) 1367 1368 Notes: 1369 This integrator can be applied to DAE. 1370 1371 Diagonally implicit general linear (DIGL) methods are a generalization of diagonally implicit Runge-Kutta (DIRK). 1372 They are represented by the tableau 1373 1374 .vb 1375 A | U 1376 ------- 1377 B | V 1378 .ve 1379 1380 combined with a vector c of abscissa. "Diagonally implicit" means that A is lower triangular. 1381 A step of the general method reads 1382 1383 .vb 1384 [ Y ] = [A U] [ Y' ] 1385 [X^k] = [B V] [X^{k-1}] 1386 .ve 1387 1388 where Y is the multivector of stage values, Y' is the multivector of stage derivatives, X^k is the Nordsieck vector of 1389 the solution at step k. The Nordsieck vector consists of the first r moments of the solution, given by 1390 1391 .vb 1392 X = [x_0,x_1,...,x_{r-1}] = [x, h x', h^2 x'', ..., h^{r-1} x^{(r-1)} ] 1393 .ve 1394 1395 If A is lower triangular, we can solve the stages (Y,Y') sequentially 1396 1397 .vb 1398 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} 1399 .ve 1400 1401 and then construct the pieces to carry to the next step 1402 1403 .vb 1404 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} 1405 .ve 1406 1407 Note that when the equations are cast in implicit form, we are using the stage equation to define y'_i 1408 in terms of y_i and known stuff (y_j for j<i and x_j for all j). 1409 1410 1411 Error estimation 1412 1413 At present, the most attractive GL methods for stiff problems are singly diagonally implicit schemes which posses 1414 Inherent Runge-Kutta Stability (IRKS). These methods have r=s, the number of items passed between steps is equal to 1415 the number of stages. The order and stage-order are one less than the number of stages. We use the error estimates 1416 in the 2007 paper which provide the following estimates 1417 1418 .vb 1419 h^{p+1} X^{(p+1)} = phi_0^T Y' + [0 psi_0^T] Xold 1420 h^{p+2} X^{(p+2)} = phi_1^T Y' + [0 psi_1^T] Xold 1421 h^{p+2} (dx'/dx) X^{(p+1)} = phi_2^T Y' + [0 psi_2^T] Xold 1422 .ve 1423 1424 These estimates are accurate to O(h^{p+3}). 1425 1426 Changing the step size 1427 1428 We use the generalized "rescale and modify" scheme, see equation (4.5) of the 2007 paper. 1429 1430 Level: beginner 1431 1432 References: 1433 + 1. - John Butcher and Z. Jackieweicz and W. Wright, On error propagation in general linear methods for 1434 ordinary differential equations, Journal of Complexity, Vol 23, 2007. 1435 - 2. - John Butcher, Numerical methods for ordinary differential equations, second edition, Wiley, 2009. 1436 1437 .seealso: TSCreate(), TS, TSSetType() 1438 1439 M*/ 1440 PETSC_EXTERN PetscErrorCode TSCreate_GLLE(TS ts) 1441 { 1442 TS_GLLE *gl; 1443 PetscErrorCode ierr; 1444 1445 PetscFunctionBegin; 1446 ierr = TSGLLEInitializePackage();CHKERRQ(ierr); 1447 1448 ierr = PetscNewLog(ts,&gl);CHKERRQ(ierr); 1449 ts->data = (void*)gl; 1450 1451 ts->ops->reset = TSReset_GLLE; 1452 ts->ops->destroy = TSDestroy_GLLE; 1453 ts->ops->view = TSView_GLLE; 1454 ts->ops->setup = TSSetUp_GLLE; 1455 ts->ops->solve = TSSolve_GLLE; 1456 ts->ops->setfromoptions = TSSetFromOptions_GLLE; 1457 ts->ops->snesfunction = SNESTSFormFunction_GLLE; 1458 ts->ops->snesjacobian = SNESTSFormJacobian_GLLE; 1459 1460 ts->usessnes = PETSC_TRUE; 1461 1462 1463 gl->max_step_rejections = 1; 1464 gl->min_order = 1; 1465 gl->max_order = 3; 1466 gl->start_order = 1; 1467 gl->current_scheme = -1; 1468 gl->extrapolate = PETSC_FALSE; 1469 1470 gl->wrms_atol = 1e-8; 1471 gl->wrms_rtol = 1e-5; 1472 1473 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSGLLESetType_C", &TSGLLESetType_GLLE);CHKERRQ(ierr); 1474 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSGLLESetAcceptType_C",&TSGLLESetAcceptType_GLLE);CHKERRQ(ierr); 1475 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSGLLEGetAdapt_C", &TSGLLEGetAdapt_GLLE);CHKERRQ(ierr); 1476 PetscFunctionReturn(0); 1477 } 1478