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