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