1 static const char help[] = "1D periodic Finite Volume solver in slope-limiter form with semidiscrete time stepping.\n" 2 "Solves scalar and vector problems, choose the physical model with -physics\n" 3 " advection - Constant coefficient scalar advection\n" 4 " u_t + (a*u)_x = 0\n" 5 " burgers - Burgers equation\n" 6 " u_t + (u^2/2)_x = 0\n" 7 " traffic - Traffic equation\n" 8 " u_t + (u*(1-u))_x = 0\n" 9 " acoustics - Acoustic wave propagation\n" 10 " u_t + (c*z*v)_x = 0\n" 11 " v_t + (c/z*u)_x = 0\n" 12 " isogas - Isothermal gas dynamics\n" 13 " rho_t + (rho*u)_x = 0\n" 14 " (rho*u)_t + (rho*u^2 + c^2*rho)_x = 0\n" 15 " shallow - Shallow water equations\n" 16 " h_t + (h*u)_x = 0\n" 17 " (h*u)_t + (h*u^2 + g*h^2/2)_x = 0\n" 18 "Some of these physical models have multiple Riemann solvers, select these with -physics_xxx_riemann\n" 19 " exact - Exact Riemann solver which usually needs to perform a Newton iteration to connect\n" 20 " the states across shocks and rarefactions\n" 21 " roe - Linearized scheme, usually with an entropy fix inside sonic rarefactions\n" 22 "The systems provide a choice of reconstructions with -physics_xxx_reconstruct\n" 23 " characteristic - Limit the characteristic variables, this is usually preferred (default)\n" 24 " conservative - Limit the conservative variables directly, can cause undesired interaction of waves\n\n" 25 "A variety of limiters for high-resolution TVD limiters are available with -limit\n" 26 " upwind,minmod,superbee,mc,vanleer,vanalbada,koren,cada-torillhon (last two are nominally third order)\n" 27 " and non-TVD schemes lax-wendroff,beam-warming,fromm\n\n" 28 "To preserve the TVD property, one should time step with a strong stability preserving method.\n" 29 "The optimal high order explicit Runge-Kutta methods in TSSSP are recommended for non-stiff problems.\n\n" 30 "Several initial conditions can be chosen with -initial N\n\n" 31 "The problem size should be set with -da_grid_x M\n\n"; 32 33 #include <petscts.h> 34 #include <petscdm.h> 35 #include <petscdmda.h> 36 #include <petscdraw.h> 37 38 #include <petsc/private/kernels/blockinvert.h> /* For the Kernel_*_gets_* stuff for BAIJ */ 39 40 PETSC_STATIC_INLINE PetscReal Sgn(PetscReal a) { return (a<0) ? -1 : 1; } 41 PETSC_STATIC_INLINE PetscReal Abs(PetscReal a) { return (a<0) ? 0 : a; } 42 PETSC_STATIC_INLINE PetscReal Sqr(PetscReal a) { return a*a; } 43 PETSC_STATIC_INLINE PetscReal MaxAbs(PetscReal a,PetscReal b) { return (PetscAbs(a) > PetscAbs(b)) ? a : b; } 44 PETSC_UNUSED PETSC_STATIC_INLINE PetscReal MinAbs(PetscReal a,PetscReal b) { return (PetscAbs(a) < PetscAbs(b)) ? a : b; } 45 PETSC_STATIC_INLINE PetscReal MinMod2(PetscReal a,PetscReal b) { return (a*b<0) ? 0 : Sgn(a)*PetscMin(PetscAbs(a),PetscAbs(b)); } 46 PETSC_STATIC_INLINE PetscReal MaxMod2(PetscReal a,PetscReal b) { return (a*b<0) ? 0 : Sgn(a)*PetscMax(PetscAbs(a),PetscAbs(b)); } 47 PETSC_STATIC_INLINE PetscReal MinMod3(PetscReal a,PetscReal b,PetscReal c) {return (a*b<0 || a*c<0) ? 0 : Sgn(a)*PetscMin(PetscAbs(a),PetscMin(PetscAbs(b),PetscAbs(c))); } 48 49 PETSC_STATIC_INLINE PetscReal RangeMod(PetscReal a,PetscReal xmin,PetscReal xmax) { PetscReal range = xmax-xmin; return xmin +PetscFmodReal(range+PetscFmodReal(a,range),range); } 50 51 /* ----------------------- Lots of limiters, these could go in a separate library ------------------------- */ 52 typedef struct _LimitInfo { 53 PetscReal hx; 54 PetscInt m; 55 } *LimitInfo; 56 static void Limit_Upwind(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt) 57 { 58 PetscInt i; 59 for (i=0; i<info->m; i++) lmt[i] = 0; 60 } 61 static void Limit_LaxWendroff(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt) 62 { 63 PetscInt i; 64 for (i=0; i<info->m; i++) lmt[i] = jR[i]; 65 } 66 static void Limit_BeamWarming(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt) 67 { 68 PetscInt i; 69 for (i=0; i<info->m; i++) lmt[i] = jL[i]; 70 } 71 static void Limit_Fromm(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt) 72 { 73 PetscInt i; 74 for (i=0; i<info->m; i++) lmt[i] = 0.5*(jL[i] + jR[i]); 75 } 76 static void Limit_Minmod(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt) 77 { 78 PetscInt i; 79 for (i=0; i<info->m; i++) lmt[i] = MinMod2(jL[i],jR[i]); 80 } 81 static void Limit_Superbee(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt) 82 { 83 PetscInt i; 84 for (i=0; i<info->m; i++) lmt[i] = MaxMod2(MinMod2(jL[i],2*jR[i]),MinMod2(2*jL[i],jR[i])); 85 } 86 static void Limit_MC(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt) 87 { 88 PetscInt i; 89 for (i=0; i<info->m; i++) lmt[i] = MinMod3(2*jL[i],0.5*(jL[i]+jR[i]),2*jR[i]); 90 } 91 static void Limit_VanLeer(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt) 92 { /* phi = (t + abs(t)) / (1 + abs(t)) */ 93 PetscInt i; 94 for (i=0; i<info->m; i++) lmt[i] = (jL[i]*Abs(jR[i]) + Abs(jL[i])*jR[i]) / (Abs(jL[i]) + Abs(jR[i]) + 1e-15); 95 } 96 static void Limit_VanAlbada(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt) /* differentiable */ 97 { /* phi = (t + t^2) / (1 + t^2) */ 98 PetscInt i; 99 for (i=0; i<info->m; i++) lmt[i] = (jL[i]*Sqr(jR[i]) + Sqr(jL[i])*jR[i]) / (Sqr(jL[i]) + Sqr(jR[i]) + 1e-15); 100 } 101 static void Limit_VanAlbadaTVD(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt) 102 { /* phi = (t + t^2) / (1 + t^2) */ 103 PetscInt i; 104 for (i=0; i<info->m; i++) lmt[i] = (jL[i]*jR[i]<0) ? 0 : (jL[i]*Sqr(jR[i]) + Sqr(jL[i])*jR[i]) / (Sqr(jL[i]) + Sqr(jR[i]) + 1e-15); 105 } 106 static void Limit_Koren(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt) /* differentiable */ 107 { /* phi = (t + 2*t^2) / (2 - t + 2*t^2) */ 108 PetscInt i; 109 for (i=0; i<info->m; i++) lmt[i] = ((jL[i]*Sqr(jR[i]) + 2*Sqr(jL[i])*jR[i])/(2*Sqr(jL[i]) - jL[i]*jR[i] + 2*Sqr(jR[i]) + 1e-15)); 110 } 111 static void Limit_KorenSym(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt) /* differentiable */ 112 { /* Symmetric version of above */ 113 PetscInt i; 114 for (i=0; i<info->m; i++) lmt[i] = (1.5*(jL[i]*Sqr(jR[i]) + Sqr(jL[i])*jR[i])/(2*Sqr(jL[i]) - jL[i]*jR[i] + 2*Sqr(jR[i]) + 1e-15)); 115 } 116 static void Limit_Koren3(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt) 117 { /* Eq 11 of Cada-Torrilhon 2009 */ 118 PetscInt i; 119 for (i=0; i<info->m; i++) lmt[i] = MinMod3(2*jL[i],(jL[i]+2*jR[i])/3,2*jR[i]); 120 } 121 static PetscReal CadaTorrilhonPhiHatR_Eq13(PetscReal L,PetscReal R) 122 { 123 return PetscMax(0,PetscMin((L+2*R)/3,PetscMax(-0.5*L,PetscMin(2*L,PetscMin((L+2*R)/3,1.6*R))))); 124 } 125 static void Limit_CadaTorrilhon2(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt) 126 { /* Cada-Torrilhon 2009, Eq 13 */ 127 PetscInt i; 128 for (i=0; i<info->m; i++) lmt[i] = CadaTorrilhonPhiHatR_Eq13(jL[i],jR[i]); 129 } 130 static void Limit_CadaTorrilhon3R(PetscReal r,LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt) 131 { /* Cada-Torrilhon 2009, Eq 22 */ 132 /* They recommend 0.001 < r < 1, but larger values are more accurate in smooth regions */ 133 const PetscReal eps = 1e-7,hx = info->hx; 134 PetscInt i; 135 for (i=0; i<info->m; i++) { 136 const PetscReal eta = (Sqr(jL[i]) + Sqr(jR[i])) / Sqr(r*hx); 137 lmt[i] = ((eta < 1-eps) ? (jL[i] + 2*jR[i]) / 3 : ((eta > 1+eps) ? CadaTorrilhonPhiHatR_Eq13(jL[i],jR[i]) : 0.5*((1-(eta-1)/eps)*(jL[i]+2*jR[i])/3 + (1+(eta+1)/eps)*CadaTorrilhonPhiHatR_Eq13(jL[i],jR[i])))); 138 } 139 } 140 static void Limit_CadaTorrilhon3R0p1(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt) 141 { 142 Limit_CadaTorrilhon3R(0.1,info,jL,jR,lmt); 143 } 144 static void Limit_CadaTorrilhon3R1(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt) 145 { 146 Limit_CadaTorrilhon3R(1,info,jL,jR,lmt); 147 } 148 static void Limit_CadaTorrilhon3R10(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt) 149 { 150 Limit_CadaTorrilhon3R(10,info,jL,jR,lmt); 151 } 152 static void Limit_CadaTorrilhon3R100(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt) 153 { 154 Limit_CadaTorrilhon3R(100,info,jL,jR,lmt); 155 } 156 157 /* --------------------------------- Finite Volume data structures ----------------------------------- */ 158 159 typedef enum {FVBC_PERIODIC, FVBC_OUTFLOW} FVBCType; 160 static const char *FVBCTypes[] = {"PERIODIC","OUTFLOW","FVBCType","FVBC_",0}; 161 typedef PetscErrorCode (*RiemannFunction)(void*,PetscInt,const PetscScalar*,const PetscScalar*,PetscScalar*,PetscReal*); 162 typedef PetscErrorCode (*ReconstructFunction)(void*,PetscInt,const PetscScalar*,PetscScalar*,PetscScalar*,PetscReal*); 163 164 typedef struct { 165 PetscErrorCode (*sample)(void*,PetscInt,FVBCType,PetscReal,PetscReal,PetscReal,PetscReal,PetscReal*); 166 RiemannFunction riemann; 167 ReconstructFunction characteristic; 168 PetscErrorCode (*destroy)(void*); 169 void *user; 170 PetscInt dof; 171 char *fieldname[16]; 172 } PhysicsCtx; 173 174 typedef struct { 175 void (*limit)(LimitInfo,const PetscScalar*,const PetscScalar*,PetscScalar*); 176 PhysicsCtx physics; 177 MPI_Comm comm; 178 char prefix[256]; 179 180 /* Local work arrays */ 181 PetscScalar *R,*Rinv; /* Characteristic basis, and it's inverse. COLUMN-MAJOR */ 182 PetscScalar *cjmpLR; /* Jumps at left and right edge of cell, in characteristic basis, len=2*dof */ 183 PetscScalar *cslope; /* Limited slope, written in characteristic basis */ 184 PetscScalar *uLR; /* Solution at left and right of interface, conservative variables, len=2*dof */ 185 PetscScalar *flux; /* Flux across interface */ 186 PetscReal *speeds; /* Speeds of each wave */ 187 188 PetscReal cfl_idt; /* Max allowable value of 1/Delta t */ 189 PetscReal cfl; 190 PetscReal xmin,xmax; 191 PetscInt initial; 192 PetscBool exact; 193 FVBCType bctype; 194 } FVCtx; 195 196 PetscErrorCode RiemannListAdd(PetscFunctionList *flist,const char *name,RiemannFunction rsolve) 197 { 198 PetscErrorCode ierr; 199 200 PetscFunctionBeginUser; 201 ierr = PetscFunctionListAdd(flist,name,rsolve);CHKERRQ(ierr); 202 PetscFunctionReturn(0); 203 } 204 205 PetscErrorCode RiemannListFind(PetscFunctionList flist,const char *name,RiemannFunction *rsolve) 206 { 207 PetscErrorCode ierr; 208 209 PetscFunctionBeginUser; 210 ierr = PetscFunctionListFind(flist,name,rsolve);CHKERRQ(ierr); 211 if (!*rsolve) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Riemann solver \"%s\" could not be found",name); 212 PetscFunctionReturn(0); 213 } 214 215 PetscErrorCode ReconstructListAdd(PetscFunctionList *flist,const char *name,ReconstructFunction r) 216 { 217 PetscErrorCode ierr; 218 219 PetscFunctionBeginUser; 220 ierr = PetscFunctionListAdd(flist,name,r);CHKERRQ(ierr); 221 PetscFunctionReturn(0); 222 } 223 224 PetscErrorCode ReconstructListFind(PetscFunctionList flist,const char *name,ReconstructFunction *r) 225 { 226 PetscErrorCode ierr; 227 228 PetscFunctionBeginUser; 229 ierr = PetscFunctionListFind(flist,name,r);CHKERRQ(ierr); 230 if (!*r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Reconstruction \"%s\" could not be found",name); 231 PetscFunctionReturn(0); 232 } 233 234 /* --------------------------------- Physics ----------------------------------- */ 235 /** 236 * Each physical model consists of Riemann solver and a function to determine the basis to use for reconstruction. These 237 * are set with the PhysicsCreate_XXX function which allocates private storage and sets these methods as well as the 238 * number of fields and their names, and a function to deallocate private storage. 239 **/ 240 241 /* First a few functions useful to several different physics */ 242 static PetscErrorCode PhysicsCharacteristic_Conservative(void *vctx,PetscInt m,const PetscScalar *u,PetscScalar *X,PetscScalar *Xi,PetscReal *speeds) 243 { 244 PetscInt i,j; 245 246 PetscFunctionBeginUser; 247 for (i=0; i<m; i++) { 248 for (j=0; j<m; j++) Xi[i*m+j] = X[i*m+j] = (PetscScalar)(i==j); 249 speeds[i] = PETSC_MAX_REAL; /* Indicates invalid */ 250 } 251 PetscFunctionReturn(0); 252 } 253 254 static PetscErrorCode PhysicsDestroy_SimpleFree(void *vctx) 255 { 256 PetscErrorCode ierr; 257 258 PetscFunctionBeginUser; 259 ierr = PetscFree(vctx);CHKERRQ(ierr); 260 PetscFunctionReturn(0); 261 } 262 263 /* --------------------------------- Advection ----------------------------------- */ 264 265 typedef struct { 266 PetscReal a; /* advective velocity */ 267 } AdvectCtx; 268 269 static PetscErrorCode PhysicsRiemann_Advect(void *vctx,PetscInt m,const PetscScalar *uL,const PetscScalar *uR,PetscScalar *flux,PetscReal *maxspeed) 270 { 271 AdvectCtx *ctx = (AdvectCtx*)vctx; 272 PetscReal speed; 273 274 PetscFunctionBeginUser; 275 speed = ctx->a; 276 flux[0] = PetscMax(0,speed)*uL[0] + PetscMin(0,speed)*uR[0]; 277 *maxspeed = speed; 278 PetscFunctionReturn(0); 279 } 280 281 static PetscErrorCode PhysicsCharacteristic_Advect(void *vctx,PetscInt m,const PetscScalar *u,PetscScalar *X,PetscScalar *Xi,PetscReal *speeds) 282 { 283 AdvectCtx *ctx = (AdvectCtx*)vctx; 284 285 PetscFunctionBeginUser; 286 X[0] = 1.; 287 Xi[0] = 1.; 288 speeds[0] = ctx->a; 289 PetscFunctionReturn(0); 290 } 291 292 static PetscErrorCode PhysicsSample_Advect(void *vctx,PetscInt initial,FVBCType bctype,PetscReal xmin,PetscReal xmax,PetscReal t,PetscReal x,PetscReal *u) 293 { 294 AdvectCtx *ctx = (AdvectCtx*)vctx; 295 PetscReal a = ctx->a,x0; 296 297 PetscFunctionBeginUser; 298 switch (bctype) { 299 case FVBC_OUTFLOW: x0 = x-a*t; break; 300 case FVBC_PERIODIC: x0 = RangeMod(x-a*t,xmin,xmax); break; 301 default: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"unknown BCType"); 302 } 303 switch (initial) { 304 case 0: u[0] = (x0 < 0) ? 1 : -1; break; 305 case 1: u[0] = (x0 < 0) ? -1 : 1; break; 306 case 2: u[0] = (0 < x0 && x0 < 1) ? 1 : 0; break; 307 case 3: u[0] = PetscSinReal(2*PETSC_PI*x0); break; 308 case 4: u[0] = PetscAbs(x0); break; 309 case 5: u[0] = (x0 < 0 || x0 > 0.5) ? 0 : PetscSqr(PetscSinReal(2*PETSC_PI*x0)); break; 310 case 6: u[0] = (x0 < 0) ? 0 : ((x0 < 1) ? x0 : ((x0 < 2) ? 2-x0 : 0)); break; 311 case 7: u[0] = PetscPowReal(PetscSinReal(PETSC_PI*x0),10.0);break; 312 default: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"unknown initial condition"); 313 } 314 PetscFunctionReturn(0); 315 } 316 317 static PetscErrorCode PhysicsCreate_Advect(FVCtx *ctx) 318 { 319 PetscErrorCode ierr; 320 AdvectCtx *user; 321 322 PetscFunctionBeginUser; 323 ierr = PetscNew(&user);CHKERRQ(ierr); 324 ctx->physics.sample = PhysicsSample_Advect; 325 ctx->physics.riemann = PhysicsRiemann_Advect; 326 ctx->physics.characteristic = PhysicsCharacteristic_Advect; 327 ctx->physics.destroy = PhysicsDestroy_SimpleFree; 328 ctx->physics.user = user; 329 ctx->physics.dof = 1; 330 ierr = PetscStrallocpy("u",&ctx->physics.fieldname[0]);CHKERRQ(ierr); 331 user->a = 1; 332 ierr = PetscOptionsBegin(ctx->comm,ctx->prefix,"Options for advection","");CHKERRQ(ierr); 333 { 334 ierr = PetscOptionsReal("-physics_advect_a","Speed","",user->a,&user->a,NULL);CHKERRQ(ierr); 335 } 336 ierr = PetscOptionsEnd();CHKERRQ(ierr); 337 PetscFunctionReturn(0); 338 } 339 340 /* --------------------------------- Burgers ----------------------------------- */ 341 342 typedef struct { 343 PetscReal lxf_speed; 344 } BurgersCtx; 345 346 static PetscErrorCode PhysicsSample_Burgers(void *vctx,PetscInt initial,FVBCType bctype,PetscReal xmin,PetscReal xmax,PetscReal t,PetscReal x,PetscReal *u) 347 { 348 PetscFunctionBeginUser; 349 if (bctype == FVBC_PERIODIC && t > 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Exact solution not implemented for periodic"); 350 switch (initial) { 351 case 0: u[0] = (x < 0) ? 1 : -1; break; 352 case 1: 353 if (x < -t) u[0] = -1; 354 else if (x < t) u[0] = x/t; 355 else u[0] = 1; 356 break; 357 case 2: 358 if (x < 0) u[0] = 0; 359 else if (x <= t) u[0] = x/t; 360 else if (x < 1+0.5*t) u[0] = 1; 361 else u[0] = 0; 362 break; 363 case 3: 364 if (x < 0.2*t) u[0] = 0.2; 365 else if (x < t) u[0] = x/t; 366 else u[0] = 1; 367 break; 368 case 4: 369 if (t > 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only initial condition available"); 370 u[0] = 0.7 + 0.3*PetscSinReal(2*PETSC_PI*((x-xmin)/(xmax-xmin))); 371 break; 372 case 5: /* Pure shock solution */ 373 if (x < 0.5*t) u[0] = 1; 374 else u[0] = 0; 375 break; 376 default: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"unknown initial condition"); 377 } 378 PetscFunctionReturn(0); 379 } 380 381 static PetscErrorCode PhysicsRiemann_Burgers_Exact(void *vctx,PetscInt m,const PetscScalar *uL,const PetscScalar *uR,PetscScalar *flux,PetscReal *maxspeed) 382 { 383 PetscFunctionBeginUser; 384 if (uL[0] < uR[0]) { /* rarefaction */ 385 flux[0] = (uL[0]*uR[0] < 0) 386 ? 0 /* sonic rarefaction */ 387 : 0.5*PetscMin(PetscSqr(uL[0]),PetscSqr(uR[0])); 388 } else { /* shock */ 389 flux[0] = 0.5*PetscMax(PetscSqr(uL[0]),PetscSqr(uR[0])); 390 } 391 *maxspeed = (PetscAbs(uL[0]) > PetscAbs(uR[0])) ? uL[0] : uR[0]; 392 PetscFunctionReturn(0); 393 } 394 395 static PetscErrorCode PhysicsRiemann_Burgers_Roe(void *vctx,PetscInt m,const PetscScalar *uL,const PetscScalar *uR,PetscScalar *flux,PetscReal *maxspeed) 396 { 397 PetscReal speed; 398 399 PetscFunctionBeginUser; 400 speed = 0.5*(uL[0] + uR[0]); 401 flux[0] = 0.25*(PetscSqr(uL[0]) + PetscSqr(uR[0])) - 0.5*PetscAbs(speed)*(uR[0]-uL[0]); 402 if (uL[0] <= 0 && 0 <= uR[0]) flux[0] = 0; /* Entropy fix for sonic rarefaction */ 403 *maxspeed = speed; 404 PetscFunctionReturn(0); 405 } 406 407 static PetscErrorCode PhysicsRiemann_Burgers_LxF(void *vctx,PetscInt m,const PetscScalar *uL,const PetscScalar *uR,PetscScalar *flux,PetscReal *maxspeed) 408 { 409 PetscReal c; 410 PetscScalar fL,fR; 411 412 PetscFunctionBeginUser; 413 c = ((BurgersCtx*)vctx)->lxf_speed; 414 fL = 0.5*PetscSqr(uL[0]); 415 fR = 0.5*PetscSqr(uR[0]); 416 flux[0] = 0.5*(fL + fR) - 0.5*c*(uR[0] - uL[0]); 417 *maxspeed = c; 418 PetscFunctionReturn(0); 419 } 420 421 static PetscErrorCode PhysicsRiemann_Burgers_Rusanov(void *vctx,PetscInt m,const PetscScalar *uL,const PetscScalar *uR,PetscScalar *flux,PetscReal *maxspeed) 422 { 423 PetscReal c; 424 PetscScalar fL,fR; 425 426 PetscFunctionBeginUser; 427 c = PetscMax(PetscAbs(uL[0]),PetscAbs(uR[0])); 428 fL = 0.5*PetscSqr(uL[0]); 429 fR = 0.5*PetscSqr(uR[0]); 430 flux[0] = 0.5*(fL + fR) - 0.5*c*(uR[0] - uL[0]); 431 *maxspeed = c; 432 PetscFunctionReturn(0); 433 } 434 435 static PetscErrorCode PhysicsCreate_Burgers(FVCtx *ctx) 436 { 437 BurgersCtx *user; 438 PetscErrorCode ierr; 439 RiemannFunction r; 440 PetscFunctionList rlist = 0; 441 char rname[256] = "exact"; 442 443 PetscFunctionBeginUser; 444 ierr = PetscNew(&user);CHKERRQ(ierr); 445 446 ctx->physics.sample = PhysicsSample_Burgers; 447 ctx->physics.characteristic = PhysicsCharacteristic_Conservative; 448 ctx->physics.destroy = PhysicsDestroy_SimpleFree; 449 ctx->physics.user = user; 450 ctx->physics.dof = 1; 451 452 ierr = PetscStrallocpy("u",&ctx->physics.fieldname[0]);CHKERRQ(ierr); 453 ierr = RiemannListAdd(&rlist,"exact", PhysicsRiemann_Burgers_Exact);CHKERRQ(ierr); 454 ierr = RiemannListAdd(&rlist,"roe", PhysicsRiemann_Burgers_Roe);CHKERRQ(ierr); 455 ierr = RiemannListAdd(&rlist,"lxf", PhysicsRiemann_Burgers_LxF);CHKERRQ(ierr); 456 ierr = RiemannListAdd(&rlist,"rusanov",PhysicsRiemann_Burgers_Rusanov);CHKERRQ(ierr); 457 ierr = PetscOptionsBegin(ctx->comm,ctx->prefix,"Options for advection","");CHKERRQ(ierr); 458 { 459 ierr = PetscOptionsFList("-physics_burgers_riemann","Riemann solver","",rlist,rname,rname,sizeof(rname),NULL);CHKERRQ(ierr); 460 } 461 ierr = PetscOptionsEnd();CHKERRQ(ierr); 462 ierr = RiemannListFind(rlist,rname,&r);CHKERRQ(ierr); 463 ierr = PetscFunctionListDestroy(&rlist);CHKERRQ(ierr); 464 ctx->physics.riemann = r; 465 466 /* * 467 * Hack to deal with LxF in semi-discrete form 468 * max speed is 1 for the basic initial conditions (where |u| <= 1) 469 * */ 470 if (r == PhysicsRiemann_Burgers_LxF) user->lxf_speed = 1; 471 PetscFunctionReturn(0); 472 } 473 474 /* --------------------------------- Traffic ----------------------------------- */ 475 476 typedef struct { 477 PetscReal lxf_speed; 478 PetscReal a; 479 } TrafficCtx; 480 481 PETSC_STATIC_INLINE PetscScalar TrafficFlux(PetscScalar a,PetscScalar u) { return a*u*(1-u); } 482 483 static PetscErrorCode PhysicsSample_Traffic(void *vctx,PetscInt initial,FVBCType bctype,PetscReal xmin,PetscReal xmax,PetscReal t,PetscReal x,PetscReal *u) 484 { 485 PetscReal a = ((TrafficCtx*)vctx)->a; 486 487 PetscFunctionBeginUser; 488 if (bctype == FVBC_PERIODIC && t > 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Exact solution not implemented for periodic"); 489 switch (initial) { 490 case 0: 491 u[0] = (-a*t < x) ? 2 : 0; break; 492 case 1: 493 if (x < PetscMin(2*a*t,0.5+a*t)) u[0] = -1; 494 else if (x < 1) u[0] = 0; 495 else u[0] = 1; 496 break; 497 case 2: 498 if (t > 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only initial condition available"); 499 u[0] = 0.7 + 0.3*PetscSinReal(2*PETSC_PI*((x-xmin)/(xmax-xmin))); 500 break; 501 default: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"unknown initial condition"); 502 } 503 PetscFunctionReturn(0); 504 } 505 506 static PetscErrorCode PhysicsRiemann_Traffic_Exact(void *vctx,PetscInt m,const PetscScalar *uL,const PetscScalar *uR,PetscScalar *flux,PetscReal *maxspeed) 507 { 508 PetscReal a = ((TrafficCtx*)vctx)->a; 509 510 PetscFunctionBeginUser; 511 if (uL[0] < uR[0]) { 512 flux[0] = PetscMin(TrafficFlux(a,uL[0]),TrafficFlux(a,uR[0])); 513 } else { 514 flux[0] = (uR[0] < 0.5 && 0.5 < uL[0]) ? TrafficFlux(a,0.5) : PetscMax(TrafficFlux(a,uL[0]),TrafficFlux(a,uR[0])); 515 } 516 *maxspeed = a*MaxAbs(1-2*uL[0],1-2*uR[0]); 517 PetscFunctionReturn(0); 518 } 519 520 static PetscErrorCode PhysicsRiemann_Traffic_Roe(void *vctx,PetscInt m,const PetscScalar *uL,const PetscScalar *uR,PetscScalar *flux,PetscReal *maxspeed) 521 { 522 PetscReal a = ((TrafficCtx*)vctx)->a; 523 PetscReal speed; 524 525 PetscFunctionBeginUser; 526 speed = a*(1 - (uL[0] + uR[0])); 527 flux[0] = 0.5*(TrafficFlux(a,uL[0]) + TrafficFlux(a,uR[0])) - 0.5*PetscAbs(speed)*(uR[0]-uL[0]); 528 *maxspeed = speed; 529 PetscFunctionReturn(0); 530 } 531 532 static PetscErrorCode PhysicsRiemann_Traffic_LxF(void *vctx,PetscInt m,const PetscScalar *uL,const PetscScalar *uR,PetscScalar *flux,PetscReal *maxspeed) 533 { 534 TrafficCtx *phys = (TrafficCtx*)vctx; 535 PetscReal a = phys->a; 536 PetscReal speed; 537 538 PetscFunctionBeginUser; 539 speed = a*(1 - (uL[0] + uR[0])); 540 flux[0] = 0.5*(TrafficFlux(a,uL[0]) + TrafficFlux(a,uR[0])) - 0.5*phys->lxf_speed*(uR[0]-uL[0]); 541 *maxspeed = speed; 542 PetscFunctionReturn(0); 543 } 544 545 static PetscErrorCode PhysicsRiemann_Traffic_Rusanov(void *vctx,PetscInt m,const PetscScalar *uL,const PetscScalar *uR,PetscScalar *flux,PetscReal *maxspeed) 546 { 547 PetscReal a = ((TrafficCtx*)vctx)->a; 548 PetscReal speed; 549 550 PetscFunctionBeginUser; 551 speed = a*PetscMax(PetscAbs(1-2*uL[0]),PetscAbs(1-2*uR[0])); 552 flux[0] = 0.5*(TrafficFlux(a,uL[0]) + TrafficFlux(a,uR[0])) - 0.5*speed*(uR[0]-uL[0]); 553 *maxspeed = speed; 554 PetscFunctionReturn(0); 555 } 556 557 static PetscErrorCode PhysicsCreate_Traffic(FVCtx *ctx) 558 { 559 PetscErrorCode ierr; 560 TrafficCtx *user; 561 RiemannFunction r; 562 PetscFunctionList rlist = 0; 563 char rname[256] = "exact"; 564 565 PetscFunctionBeginUser; 566 ierr = PetscNew(&user);CHKERRQ(ierr); 567 ctx->physics.sample = PhysicsSample_Traffic; 568 ctx->physics.characteristic = PhysicsCharacteristic_Conservative; 569 ctx->physics.destroy = PhysicsDestroy_SimpleFree; 570 ctx->physics.user = user; 571 ctx->physics.dof = 1; 572 573 ierr = PetscStrallocpy("density",&ctx->physics.fieldname[0]);CHKERRQ(ierr); 574 user->a = 0.5; 575 ierr = RiemannListAdd(&rlist,"exact", PhysicsRiemann_Traffic_Exact);CHKERRQ(ierr); 576 ierr = RiemannListAdd(&rlist,"roe", PhysicsRiemann_Traffic_Roe);CHKERRQ(ierr); 577 ierr = RiemannListAdd(&rlist,"lxf", PhysicsRiemann_Traffic_LxF);CHKERRQ(ierr); 578 ierr = RiemannListAdd(&rlist,"rusanov",PhysicsRiemann_Traffic_Rusanov);CHKERRQ(ierr); 579 ierr = PetscOptionsBegin(ctx->comm,ctx->prefix,"Options for Traffic","");CHKERRQ(ierr); 580 ierr = PetscOptionsReal("-physics_traffic_a","Flux = a*u*(1-u)","",user->a,&user->a,NULL);CHKERRQ(ierr); 581 ierr = PetscOptionsFList("-physics_traffic_riemann","Riemann solver","",rlist,rname,rname,sizeof(rname),NULL);CHKERRQ(ierr); 582 ierr = PetscOptionsEnd();CHKERRQ(ierr); 583 584 ierr = RiemannListFind(rlist,rname,&r);CHKERRQ(ierr); 585 ierr = PetscFunctionListDestroy(&rlist);CHKERRQ(ierr); 586 587 ctx->physics.riemann = r; 588 589 /* * 590 * Hack to deal with LxF in semi-discrete form 591 * max speed is 3*a for the basic initial conditions (-1 <= u <= 2) 592 * */ 593 if (r == PhysicsRiemann_Traffic_LxF) user->lxf_speed = 3*user->a; 594 PetscFunctionReturn(0); 595 } 596 597 /* --------------------------------- Linear Acoustics ----------------------------------- */ 598 599 /* Flux: u_t + (A u)_x 600 * z = sqrt(rho*bulk), c = sqrt(rho/bulk) 601 * Spectral decomposition: A = R * D * Rinv 602 * [ cz] = [-z z] [-c ] [-1/2z 1/2] 603 * [c/z ] = [ 1 1] [ c] [ 1/2z 1/2] 604 * 605 * We decompose this into the left-traveling waves Al = R * D^- Rinv 606 * and the right-traveling waves Ar = R * D^+ * Rinv 607 * Multiplying out these expressions produces the following two matrices 608 */ 609 610 typedef struct { 611 PetscReal c; /* speed of sound: c = sqrt(bulk/rho) */ 612 PetscReal z; /* impedence: z = sqrt(rho*bulk) */ 613 } AcousticsCtx; 614 615 PETSC_UNUSED PETSC_STATIC_INLINE void AcousticsFlux(AcousticsCtx *ctx,const PetscScalar *u,PetscScalar *f) 616 { 617 f[0] = ctx->c*ctx->z*u[1]; 618 f[1] = ctx->c/ctx->z*u[0]; 619 } 620 621 static PetscErrorCode PhysicsCharacteristic_Acoustics(void *vctx,PetscInt m,const PetscScalar *u,PetscScalar *X,PetscScalar *Xi,PetscReal *speeds) 622 { 623 AcousticsCtx *phys = (AcousticsCtx*)vctx; 624 PetscReal z = phys->z,c = phys->c; 625 626 PetscFunctionBeginUser; 627 X[0*2+0] = -z; 628 X[0*2+1] = z; 629 X[1*2+0] = 1; 630 X[1*2+1] = 1; 631 Xi[0*2+0] = -1./(2*z); 632 Xi[0*2+1] = 1./2; 633 Xi[1*2+0] = 1./(2*z); 634 Xi[1*2+1] = 1./2; 635 speeds[0] = -c; 636 speeds[1] = c; 637 PetscFunctionReturn(0); 638 } 639 640 static PetscErrorCode PhysicsSample_Acoustics_Initial(AcousticsCtx *phys,PetscInt initial,PetscReal xmin,PetscReal xmax,PetscReal x,PetscReal *u) 641 { 642 PetscFunctionBeginUser; 643 switch (initial) { 644 case 0: 645 u[0] = (PetscAbs((x - xmin)/(xmax - xmin) - 0.2) < 0.1) ? 1 : 0.5; 646 u[1] = (PetscAbs((x - xmin)/(xmax - xmin) - 0.7) < 0.1) ? 1 : -0.5; 647 break; 648 case 1: 649 u[0] = PetscCosReal(3 * 2*PETSC_PI*x/(xmax-xmin)); 650 u[1] = PetscExpReal(-PetscSqr(x - (xmax + xmin)/2) / (2*PetscSqr(0.2*(xmax - xmin)))) - 0.5; 651 break; 652 default: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"unknown initial condition"); 653 } 654 PetscFunctionReturn(0); 655 } 656 657 static PetscErrorCode PhysicsSample_Acoustics(void *vctx,PetscInt initial,FVBCType bctype,PetscReal xmin,PetscReal xmax,PetscReal t,PetscReal x,PetscReal *u) 658 { 659 AcousticsCtx *phys = (AcousticsCtx*)vctx; 660 PetscReal c = phys->c; 661 PetscReal x0a,x0b,u0a[2],u0b[2],tmp[2]; 662 PetscReal X[2][2],Xi[2][2],dummy[2]; 663 PetscErrorCode ierr; 664 665 PetscFunctionBeginUser; 666 switch (bctype) { 667 case FVBC_OUTFLOW: 668 x0a = x+c*t; 669 x0b = x-c*t; 670 break; 671 case FVBC_PERIODIC: 672 x0a = RangeMod(x+c*t,xmin,xmax); 673 x0b = RangeMod(x-c*t,xmin,xmax); 674 break; 675 default: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"unknown BCType"); 676 } 677 ierr = PhysicsSample_Acoustics_Initial(phys,initial,xmin,xmax,x0a,u0a);CHKERRQ(ierr); 678 ierr = PhysicsSample_Acoustics_Initial(phys,initial,xmin,xmax,x0b,u0b);CHKERRQ(ierr); 679 ierr = PhysicsCharacteristic_Acoustics(vctx,2,u,&X[0][0],&Xi[0][0],dummy);CHKERRQ(ierr); 680 tmp[0] = Xi[0][0]*u0a[0] + Xi[0][1]*u0a[1]; 681 tmp[1] = Xi[1][0]*u0b[0] + Xi[1][1]*u0b[1]; 682 u[0] = X[0][0]*tmp[0] + X[0][1]*tmp[1]; 683 u[1] = X[1][0]*tmp[0] + X[1][1]*tmp[1]; 684 PetscFunctionReturn(0); 685 } 686 687 static PetscErrorCode PhysicsRiemann_Acoustics_Exact(void *vctx,PetscInt m,const PetscScalar *uL,const PetscScalar *uR,PetscScalar *flux,PetscReal *maxspeed) 688 { 689 AcousticsCtx *phys = (AcousticsCtx*)vctx; 690 PetscReal c = phys->c,z = phys->z; 691 PetscReal 692 Al[2][2] = {{-c/2 , c*z/2 }, 693 {c/(2*z) , -c/2 }}, /* Left traveling waves */ 694 Ar[2][2] = {{c/2 , c*z/2 }, 695 {c/(2*z) , c/2 }}; /* Right traveling waves */ 696 697 PetscFunctionBeginUser; 698 flux[0] = Al[0][0]*uR[0] + Al[0][1]*uR[1] + Ar[0][0]*uL[0] + Ar[0][1]*uL[1]; 699 flux[1] = Al[1][0]*uR[0] + Al[1][1]*uR[1] + Ar[1][0]*uL[0] + Ar[1][1]*uL[1]; 700 *maxspeed = c; 701 PetscFunctionReturn(0); 702 } 703 704 static PetscErrorCode PhysicsCreate_Acoustics(FVCtx *ctx) 705 { 706 PetscErrorCode ierr; 707 AcousticsCtx *user; 708 PetscFunctionList rlist = 0,rclist = 0; 709 char rname[256] = "exact",rcname[256] = "characteristic"; 710 711 PetscFunctionBeginUser; 712 ierr = PetscNew(&user);CHKERRQ(ierr); 713 ctx->physics.sample = PhysicsSample_Acoustics; 714 ctx->physics.destroy = PhysicsDestroy_SimpleFree; 715 ctx->physics.user = user; 716 ctx->physics.dof = 2; 717 718 ierr = PetscStrallocpy("u",&ctx->physics.fieldname[0]);CHKERRQ(ierr); 719 ierr = PetscStrallocpy("v",&ctx->physics.fieldname[1]);CHKERRQ(ierr); 720 721 user->c = 1; 722 user->z = 1; 723 724 ierr = RiemannListAdd(&rlist,"exact", PhysicsRiemann_Acoustics_Exact);CHKERRQ(ierr); 725 ierr = ReconstructListAdd(&rclist,"characteristic",PhysicsCharacteristic_Acoustics);CHKERRQ(ierr); 726 ierr = ReconstructListAdd(&rclist,"conservative",PhysicsCharacteristic_Conservative);CHKERRQ(ierr); 727 ierr = PetscOptionsBegin(ctx->comm,ctx->prefix,"Options for linear Acoustics","");CHKERRQ(ierr); 728 { 729 ierr = PetscOptionsReal("-physics_acoustics_c","c = sqrt(bulk/rho)","",user->c,&user->c,NULL);CHKERRQ(ierr); 730 ierr = PetscOptionsReal("-physics_acoustics_z","z = sqrt(bulk*rho)","",user->z,&user->z,NULL);CHKERRQ(ierr); 731 ierr = PetscOptionsFList("-physics_acoustics_riemann","Riemann solver","",rlist,rname,rname,sizeof(rname),NULL);CHKERRQ(ierr); 732 ierr = PetscOptionsFList("-physics_acoustics_reconstruct","Reconstruction","",rclist,rcname,rcname,sizeof(rcname),NULL);CHKERRQ(ierr); 733 } 734 ierr = PetscOptionsEnd();CHKERRQ(ierr); 735 ierr = RiemannListFind(rlist,rname,&ctx->physics.riemann);CHKERRQ(ierr); 736 ierr = ReconstructListFind(rclist,rcname,&ctx->physics.characteristic);CHKERRQ(ierr); 737 ierr = PetscFunctionListDestroy(&rlist);CHKERRQ(ierr); 738 ierr = PetscFunctionListDestroy(&rclist);CHKERRQ(ierr); 739 PetscFunctionReturn(0); 740 } 741 742 /* --------------------------------- Isothermal Gas Dynamics ----------------------------------- */ 743 744 typedef struct { 745 PetscReal acoustic_speed; 746 } IsoGasCtx; 747 748 PETSC_STATIC_INLINE void IsoGasFlux(PetscReal c,const PetscScalar *u,PetscScalar *f) 749 { 750 f[0] = u[1]; 751 f[1] = PetscSqr(u[1])/u[0] + c*c*u[0]; 752 } 753 754 static PetscErrorCode PhysicsSample_IsoGas(void *vctx,PetscInt initial,FVBCType bctype,PetscReal xmin,PetscReal xmax,PetscReal t,PetscReal x,PetscReal *u) 755 { 756 PetscFunctionBeginUser; 757 if (t > 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Exact solutions not implemented for t > 0"); 758 switch (initial) { 759 case 0: 760 u[0] = (x < 0) ? 1 : 0.5; 761 u[1] = (x < 0) ? 1 : 0.7; 762 break; 763 case 1: 764 u[0] = 1+0.5*PetscSinReal(2*PETSC_PI*x); 765 u[1] = 1*u[0]; 766 break; 767 default: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"unknown initial condition"); 768 } 769 PetscFunctionReturn(0); 770 } 771 772 static PetscErrorCode PhysicsRiemann_IsoGas_Roe(void *vctx,PetscInt m,const PetscScalar *uL,const PetscScalar *uR,PetscScalar *flux,PetscReal *maxspeed) 773 { 774 IsoGasCtx *phys = (IsoGasCtx*)vctx; 775 PetscReal c = phys->acoustic_speed; 776 PetscScalar ubar,du[2],a[2],fL[2],fR[2],lam[2],ustar[2],R[2][2]; 777 PetscInt i; 778 779 PetscFunctionBeginUser; 780 ubar = (uL[1]/PetscSqrtScalar(uL[0]) + uR[1]/PetscSqrtScalar(uR[0])) / (PetscSqrtScalar(uL[0]) + PetscSqrtScalar(uR[0])); 781 /* write fluxuations in characteristic basis */ 782 du[0] = uR[0] - uL[0]; 783 du[1] = uR[1] - uL[1]; 784 a[0] = (1/(2*c)) * ((ubar + c)*du[0] - du[1]); 785 a[1] = (1/(2*c)) * ((-ubar + c)*du[0] + du[1]); 786 /* wave speeds */ 787 lam[0] = ubar - c; 788 lam[1] = ubar + c; 789 /* Right eigenvectors */ 790 R[0][0] = 1; R[0][1] = ubar-c; 791 R[1][0] = 1; R[1][1] = ubar+c; 792 /* Compute state in star region (between the 1-wave and 2-wave) */ 793 for (i=0; i<2; i++) ustar[i] = uL[i] + a[0]*R[0][i]; 794 if (uL[1]/uL[0] < c && c < ustar[1]/ustar[0]) { /* 1-wave is sonic rarefaction */ 795 PetscScalar ufan[2]; 796 ufan[0] = uL[0]*PetscExpScalar(uL[1]/(uL[0]*c) - 1); 797 ufan[1] = c*ufan[0]; 798 IsoGasFlux(c,ufan,flux); 799 } else if (ustar[1]/ustar[0] < -c && -c < uR[1]/uR[0]) { /* 2-wave is sonic rarefaction */ 800 PetscScalar ufan[2]; 801 ufan[0] = uR[0]*PetscExpScalar(-uR[1]/(uR[0]*c) - 1); 802 ufan[1] = -c*ufan[0]; 803 IsoGasFlux(c,ufan,flux); 804 } else { /* Centered form */ 805 IsoGasFlux(c,uL,fL); 806 IsoGasFlux(c,uR,fR); 807 for (i=0; i<2; i++) { 808 PetscScalar absdu = PetscAbsScalar(lam[0])*a[0]*R[0][i] + PetscAbsScalar(lam[1])*a[1]*R[1][i]; 809 flux[i] = 0.5*(fL[i]+fR[i]) - 0.5*absdu; 810 } 811 } 812 *maxspeed = MaxAbs(lam[0],lam[1]); 813 PetscFunctionReturn(0); 814 } 815 816 static PetscErrorCode PhysicsRiemann_IsoGas_Exact(void *vctx,PetscInt m,const PetscScalar *uL,const PetscScalar *uR,PetscScalar *flux,PetscReal *maxspeed) 817 { 818 IsoGasCtx *phys = (IsoGasCtx*)vctx; 819 PetscReal c = phys->acoustic_speed; 820 PetscScalar ustar[2]; 821 struct {PetscScalar rho,u;} L = {uL[0],uL[1]/uL[0]},R = {uR[0],uR[1]/uR[0]},star; 822 PetscInt i; 823 824 PetscFunctionBeginUser; 825 if (!(L.rho > 0 && R.rho > 0)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Reconstructed density is negative"); 826 { 827 /* Solve for star state */ 828 PetscScalar res,tmp,rho = 0.5*(L.rho + R.rho); /* initial guess */ 829 for (i=0; i<20; i++) { 830 PetscScalar fr,fl,dfr,dfl; 831 fl = (L.rho < rho) 832 ? (rho-L.rho)/PetscSqrtScalar(L.rho*rho) /* shock */ 833 : PetscLogScalar(rho) - PetscLogScalar(L.rho); /* rarefaction */ 834 fr = (R.rho < rho) 835 ? (rho-R.rho)/PetscSqrtScalar(R.rho*rho) /* shock */ 836 : PetscLogScalar(rho) - PetscLogScalar(R.rho); /* rarefaction */ 837 res = R.u-L.u + c*(fr+fl); 838 if (PetscIsInfOrNanScalar(res)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"Infinity or Not-a-Number generated in computation"); 839 if (PetscAbsScalar(res) < 1e-10) { 840 star.rho = rho; 841 star.u = L.u - c*fl; 842 goto converged; 843 } 844 dfl = (L.rho < rho) ? 1/PetscSqrtScalar(L.rho*rho)*(1 - 0.5*(rho-L.rho)/rho) : 1/rho; 845 dfr = (R.rho < rho) ? 1/PetscSqrtScalar(R.rho*rho)*(1 - 0.5*(rho-R.rho)/rho) : 1/rho; 846 tmp = rho - res/(c*(dfr+dfl)); 847 if (tmp <= 0) rho /= 2; /* Guard against Newton shooting off to a negative density */ 848 else rho = tmp; 849 if (!((rho > 0) && PetscIsNormalScalar(rho))) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_FP,"non-normal iterate rho=%g",(double)PetscRealPart(rho)); 850 } 851 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_CONV_FAILED,"Newton iteration for star.rho diverged after %D iterations",i); 852 } 853 converged: 854 if (L.u-c < 0 && 0 < star.u-c) { /* 1-wave is sonic rarefaction */ 855 PetscScalar ufan[2]; 856 ufan[0] = L.rho*PetscExpScalar(L.u/c - 1); 857 ufan[1] = c*ufan[0]; 858 IsoGasFlux(c,ufan,flux); 859 } else if (star.u+c < 0 && 0 < R.u+c) { /* 2-wave is sonic rarefaction */ 860 PetscScalar ufan[2]; 861 ufan[0] = R.rho*PetscExpScalar(-R.u/c - 1); 862 ufan[1] = -c*ufan[0]; 863 IsoGasFlux(c,ufan,flux); 864 } else if ((L.rho >= star.rho && L.u-c >= 0) || (L.rho < star.rho && (star.rho*star.u-L.rho*L.u)/(star.rho-L.rho) > 0)) { 865 /* 1-wave is supersonic rarefaction, or supersonic shock */ 866 IsoGasFlux(c,uL,flux); 867 } else if ((star.rho <= R.rho && R.u+c <= 0) || (star.rho > R.rho && (R.rho*R.u-star.rho*star.u)/(R.rho-star.rho) < 0)) { 868 /* 2-wave is supersonic rarefaction or supersonic shock */ 869 IsoGasFlux(c,uR,flux); 870 } else { 871 ustar[0] = star.rho; 872 ustar[1] = star.rho*star.u; 873 IsoGasFlux(c,ustar,flux); 874 } 875 *maxspeed = MaxAbs(MaxAbs(star.u-c,star.u+c),MaxAbs(L.u-c,R.u+c)); 876 PetscFunctionReturn(0); 877 } 878 879 static PetscErrorCode PhysicsRiemann_IsoGas_Rusanov(void *vctx,PetscInt m,const PetscScalar *uL,const PetscScalar *uR,PetscScalar *flux,PetscReal *maxspeed) 880 { 881 IsoGasCtx *phys = (IsoGasCtx*)vctx; 882 PetscScalar c = phys->acoustic_speed,fL[2],fR[2],s; 883 struct {PetscScalar rho,u;} L = {uL[0],uL[1]/uL[0]},R = {uR[0],uR[1]/uR[0]}; 884 885 PetscFunctionBeginUser; 886 if (!(L.rho > 0 && R.rho > 0)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Reconstructed density is negative"); 887 IsoGasFlux(c,uL,fL); 888 IsoGasFlux(c,uR,fR); 889 s = PetscMax(PetscAbs(L.u),PetscAbs(R.u))+c; 890 flux[0] = 0.5*(fL[0] + fR[0]) + 0.5*s*(uL[0] - uR[0]); 891 flux[1] = 0.5*(fL[1] + fR[1]) + 0.5*s*(uL[1] - uR[1]); 892 *maxspeed = s; 893 PetscFunctionReturn(0); 894 } 895 896 static PetscErrorCode PhysicsCharacteristic_IsoGas(void *vctx,PetscInt m,const PetscScalar *u,PetscScalar *X,PetscScalar *Xi,PetscReal *speeds) 897 { 898 IsoGasCtx *phys = (IsoGasCtx*)vctx; 899 PetscReal c = phys->acoustic_speed; 900 PetscErrorCode ierr; 901 902 PetscFunctionBeginUser; 903 speeds[0] = u[1]/u[0] - c; 904 speeds[1] = u[1]/u[0] + c; 905 X[0*2+0] = 1; 906 X[0*2+1] = speeds[0]; 907 X[1*2+0] = 1; 908 X[1*2+1] = speeds[1]; 909 ierr = PetscArraycpy(Xi,X,4);CHKERRQ(ierr); 910 ierr = PetscKernel_A_gets_inverse_A_2(Xi,0,PETSC_FALSE,NULL);CHKERRQ(ierr); 911 PetscFunctionReturn(0); 912 } 913 914 static PetscErrorCode PhysicsCreate_IsoGas(FVCtx *ctx) 915 { 916 PetscErrorCode ierr; 917 IsoGasCtx *user; 918 PetscFunctionList rlist = 0,rclist = 0; 919 char rname[256] = "exact",rcname[256] = "characteristic"; 920 921 PetscFunctionBeginUser; 922 ierr = PetscNew(&user);CHKERRQ(ierr); 923 ctx->physics.sample = PhysicsSample_IsoGas; 924 ctx->physics.destroy = PhysicsDestroy_SimpleFree; 925 ctx->physics.user = user; 926 ctx->physics.dof = 2; 927 928 ierr = PetscStrallocpy("density",&ctx->physics.fieldname[0]);CHKERRQ(ierr); 929 ierr = PetscStrallocpy("momentum",&ctx->physics.fieldname[1]);CHKERRQ(ierr); 930 931 user->acoustic_speed = 1; 932 933 ierr = RiemannListAdd(&rlist,"exact", PhysicsRiemann_IsoGas_Exact);CHKERRQ(ierr); 934 ierr = RiemannListAdd(&rlist,"roe", PhysicsRiemann_IsoGas_Roe);CHKERRQ(ierr); 935 ierr = RiemannListAdd(&rlist,"rusanov",PhysicsRiemann_IsoGas_Rusanov);CHKERRQ(ierr); 936 ierr = ReconstructListAdd(&rclist,"characteristic",PhysicsCharacteristic_IsoGas);CHKERRQ(ierr); 937 ierr = ReconstructListAdd(&rclist,"conservative",PhysicsCharacteristic_Conservative);CHKERRQ(ierr); 938 ierr = PetscOptionsBegin(ctx->comm,ctx->prefix,"Options for IsoGas","");CHKERRQ(ierr); 939 ierr = PetscOptionsReal("-physics_isogas_acoustic_speed","Acoustic speed","",user->acoustic_speed,&user->acoustic_speed,NULL);CHKERRQ(ierr); 940 ierr = PetscOptionsFList("-physics_isogas_riemann","Riemann solver","",rlist,rname,rname,sizeof(rname),NULL);CHKERRQ(ierr); 941 ierr = PetscOptionsFList("-physics_isogas_reconstruct","Reconstruction","",rclist,rcname,rcname,sizeof(rcname),NULL);CHKERRQ(ierr); 942 ierr = PetscOptionsEnd();CHKERRQ(ierr); 943 ierr = RiemannListFind(rlist,rname,&ctx->physics.riemann);CHKERRQ(ierr); 944 ierr = ReconstructListFind(rclist,rcname,&ctx->physics.characteristic);CHKERRQ(ierr); 945 ierr = PetscFunctionListDestroy(&rlist);CHKERRQ(ierr); 946 ierr = PetscFunctionListDestroy(&rclist);CHKERRQ(ierr); 947 PetscFunctionReturn(0); 948 } 949 950 /* --------------------------------- Shallow Water ----------------------------------- */ 951 typedef struct { 952 PetscReal gravity; 953 } ShallowCtx; 954 955 PETSC_STATIC_INLINE void ShallowFlux(ShallowCtx *phys,const PetscScalar *u,PetscScalar *f) 956 { 957 f[0] = u[1]; 958 f[1] = PetscSqr(u[1])/u[0] + 0.5*phys->gravity*PetscSqr(u[0]); 959 } 960 961 static PetscErrorCode PhysicsRiemann_Shallow_Exact(void *vctx,PetscInt m,const PetscScalar *uL,const PetscScalar *uR,PetscScalar *flux,PetscReal *maxspeed) 962 { 963 ShallowCtx *phys = (ShallowCtx*)vctx; 964 PetscScalar g = phys->gravity,ustar[2],cL,cR,c,cstar; 965 struct {PetscScalar h,u;} L = {uL[0],uL[1]/uL[0]},R = {uR[0],uR[1]/uR[0]},star; 966 PetscInt i; 967 968 PetscFunctionBeginUser; 969 if (!(L.h > 0 && R.h > 0)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Reconstructed thickness is negative"); 970 cL = PetscSqrtScalar(g*L.h); 971 cR = PetscSqrtScalar(g*R.h); 972 c = PetscMax(cL,cR); 973 { 974 /* Solve for star state */ 975 const PetscInt maxits = 50; 976 PetscScalar tmp,res,res0=0,h0,h = 0.5*(L.h + R.h); /* initial guess */ 977 h0 = h; 978 for (i=0; i<maxits; i++) { 979 PetscScalar fr,fl,dfr,dfl; 980 fl = (L.h < h) 981 ? PetscSqrtScalar(0.5*g*(h*h - L.h*L.h)*(1/L.h - 1/h)) /* shock */ 982 : 2*PetscSqrtScalar(g*h) - 2*PetscSqrtScalar(g*L.h); /* rarefaction */ 983 fr = (R.h < h) 984 ? PetscSqrtScalar(0.5*g*(h*h - R.h*R.h)*(1/R.h - 1/h)) /* shock */ 985 : 2*PetscSqrtScalar(g*h) - 2*PetscSqrtScalar(g*R.h); /* rarefaction */ 986 res = R.u - L.u + fr + fl; 987 if (PetscIsInfOrNanScalar(res)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"Infinity or Not-a-Number generated in computation"); 988 if (PetscAbsScalar(res) < 1e-8 || (i > 0 && PetscAbsScalar(h-h0) < 1e-8)) { 989 star.h = h; 990 star.u = L.u - fl; 991 goto converged; 992 } else if (i > 0 && PetscAbsScalar(res) >= PetscAbsScalar(res0)) { /* Line search */ 993 h = 0.8*h0 + 0.2*h; 994 continue; 995 } 996 /* Accept the last step and take another */ 997 res0 = res; 998 h0 = h; 999 dfl = (L.h < h) ? 0.5/fl*0.5*g*(-L.h*L.h/(h*h) - 1 + 2*h/L.h) : PetscSqrtScalar(g/h); 1000 dfr = (R.h < h) ? 0.5/fr*0.5*g*(-R.h*R.h/(h*h) - 1 + 2*h/R.h) : PetscSqrtScalar(g/h); 1001 tmp = h - res/(dfr+dfl); 1002 if (tmp <= 0) h /= 2; /* Guard against Newton shooting off to a negative thickness */ 1003 else h = tmp; 1004 if (!((h > 0) && PetscIsNormalScalar(h))) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_FP,"non-normal iterate h=%g",(double)h); 1005 } 1006 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_CONV_FAILED,"Newton iteration for star.h diverged after %D iterations",i); 1007 } 1008 converged: 1009 cstar = PetscSqrtScalar(g*star.h); 1010 if (L.u-cL < 0 && 0 < star.u-cstar) { /* 1-wave is sonic rarefaction */ 1011 PetscScalar ufan[2]; 1012 ufan[0] = 1/g*PetscSqr(L.u/3 + 2./3*cL); 1013 ufan[1] = PetscSqrtScalar(g*ufan[0])*ufan[0]; 1014 ShallowFlux(phys,ufan,flux); 1015 } else if (star.u+cstar < 0 && 0 < R.u+cR) { /* 2-wave is sonic rarefaction */ 1016 PetscScalar ufan[2]; 1017 ufan[0] = 1/g*PetscSqr(R.u/3 - 2./3*cR); 1018 ufan[1] = -PetscSqrtScalar(g*ufan[0])*ufan[0]; 1019 ShallowFlux(phys,ufan,flux); 1020 } else if ((L.h >= star.h && L.u-c >= 0) || (L.h<star.h && (star.h*star.u-L.h*L.u)/(star.h-L.h) > 0)) { 1021 /* 1-wave is right-travelling shock (supersonic) */ 1022 ShallowFlux(phys,uL,flux); 1023 } else if ((star.h <= R.h && R.u+c <= 0) || (star.h>R.h && (R.h*R.u-star.h*star.h)/(R.h-star.h) < 0)) { 1024 /* 2-wave is left-travelling shock (supersonic) */ 1025 ShallowFlux(phys,uR,flux); 1026 } else { 1027 ustar[0] = star.h; 1028 ustar[1] = star.h*star.u; 1029 ShallowFlux(phys,ustar,flux); 1030 } 1031 *maxspeed = MaxAbs(MaxAbs(star.u-cstar,star.u+cstar),MaxAbs(L.u-cL,R.u+cR)); 1032 PetscFunctionReturn(0); 1033 } 1034 1035 static PetscErrorCode PhysicsRiemann_Shallow_Rusanov(void *vctx,PetscInt m,const PetscScalar *uL,const PetscScalar *uR,PetscScalar *flux,PetscReal *maxspeed) 1036 { 1037 ShallowCtx *phys = (ShallowCtx*)vctx; 1038 PetscScalar g = phys->gravity,fL[2],fR[2],s; 1039 struct {PetscScalar h,u;} L = {uL[0],uL[1]/uL[0]},R = {uR[0],uR[1]/uR[0]}; 1040 1041 PetscFunctionBeginUser; 1042 if (!(L.h > 0 && R.h > 0)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Reconstructed thickness is negative"); 1043 ShallowFlux(phys,uL,fL); 1044 ShallowFlux(phys,uR,fR); 1045 s = PetscMax(PetscAbs(L.u)+PetscSqrtScalar(g*L.h),PetscAbs(R.u)+PetscSqrtScalar(g*R.h)); 1046 flux[0] = 0.5*(fL[0] + fR[0]) + 0.5*s*(uL[0] - uR[0]); 1047 flux[1] = 0.5*(fL[1] + fR[1]) + 0.5*s*(uL[1] - uR[1]); 1048 *maxspeed = s; 1049 PetscFunctionReturn(0); 1050 } 1051 1052 static PetscErrorCode PhysicsCharacteristic_Shallow(void *vctx,PetscInt m,const PetscScalar *u,PetscScalar *X,PetscScalar *Xi,PetscReal *speeds) 1053 { 1054 ShallowCtx *phys = (ShallowCtx*)vctx; 1055 PetscReal c; 1056 PetscErrorCode ierr; 1057 1058 PetscFunctionBeginUser; 1059 c = PetscSqrtScalar(u[0]*phys->gravity); 1060 speeds[0] = u[1]/u[0] - c; 1061 speeds[1] = u[1]/u[0] + c; 1062 X[0*2+0] = 1; 1063 X[0*2+1] = speeds[0]; 1064 X[1*2+0] = 1; 1065 X[1*2+1] = speeds[1]; 1066 ierr = PetscArraycpy(Xi,X,4);CHKERRQ(ierr); 1067 ierr = PetscKernel_A_gets_inverse_A_2(Xi,0,PETSC_FALSE,NULL);CHKERRQ(ierr); 1068 PetscFunctionReturn(0); 1069 } 1070 1071 static PetscErrorCode PhysicsCreate_Shallow(FVCtx *ctx) 1072 { 1073 PetscErrorCode ierr; 1074 ShallowCtx *user; 1075 PetscFunctionList rlist = 0,rclist = 0; 1076 char rname[256] = "exact",rcname[256] = "characteristic"; 1077 1078 PetscFunctionBeginUser; 1079 ierr = PetscNew(&user);CHKERRQ(ierr); 1080 /* Shallow water and Isothermal Gas dynamics are similar so we reuse initial conditions for now */ 1081 ctx->physics.sample = PhysicsSample_IsoGas; 1082 ctx->physics.destroy = PhysicsDestroy_SimpleFree; 1083 ctx->physics.user = user; 1084 ctx->physics.dof = 2; 1085 1086 ierr = PetscStrallocpy("density",&ctx->physics.fieldname[0]);CHKERRQ(ierr); 1087 ierr = PetscStrallocpy("momentum",&ctx->physics.fieldname[1]);CHKERRQ(ierr); 1088 1089 user->gravity = 1; 1090 1091 ierr = RiemannListAdd(&rlist,"exact", PhysicsRiemann_Shallow_Exact);CHKERRQ(ierr); 1092 ierr = RiemannListAdd(&rlist,"rusanov",PhysicsRiemann_Shallow_Rusanov);CHKERRQ(ierr); 1093 ierr = ReconstructListAdd(&rclist,"characteristic",PhysicsCharacteristic_Shallow);CHKERRQ(ierr); 1094 ierr = ReconstructListAdd(&rclist,"conservative",PhysicsCharacteristic_Conservative);CHKERRQ(ierr); 1095 ierr = PetscOptionsBegin(ctx->comm,ctx->prefix,"Options for Shallow","");CHKERRQ(ierr); 1096 ierr = PetscOptionsReal("-physics_shallow_gravity","Gravity","",user->gravity,&user->gravity,NULL);CHKERRQ(ierr); 1097 ierr = PetscOptionsFList("-physics_shallow_riemann","Riemann solver","",rlist,rname,rname,sizeof(rname),NULL);CHKERRQ(ierr); 1098 ierr = PetscOptionsFList("-physics_shallow_reconstruct","Reconstruction","",rclist,rcname,rcname,sizeof(rcname),NULL);CHKERRQ(ierr); 1099 ierr = PetscOptionsEnd();CHKERRQ(ierr); 1100 ierr = RiemannListFind(rlist,rname,&ctx->physics.riemann);CHKERRQ(ierr); 1101 ierr = ReconstructListFind(rclist,rcname,&ctx->physics.characteristic);CHKERRQ(ierr); 1102 ierr = PetscFunctionListDestroy(&rlist);CHKERRQ(ierr); 1103 ierr = PetscFunctionListDestroy(&rclist);CHKERRQ(ierr); 1104 PetscFunctionReturn(0); 1105 } 1106 1107 /* --------------------------------- Finite Volume Solver ----------------------------------- */ 1108 1109 static PetscErrorCode FVRHSFunction(TS ts,PetscReal time,Vec X,Vec F,void *vctx) 1110 { 1111 FVCtx *ctx = (FVCtx*)vctx; 1112 PetscErrorCode ierr; 1113 PetscInt i,j,k,Mx,dof,xs,xm; 1114 PetscReal hx,cfl_idt = 0; 1115 PetscScalar *x,*f,*slope; 1116 Vec Xloc; 1117 DM da; 1118 1119 PetscFunctionBeginUser; 1120 ierr = TSGetDM(ts,&da);CHKERRQ(ierr); 1121 ierr = DMGetLocalVector(da,&Xloc);CHKERRQ(ierr); 1122 ierr = DMDAGetInfo(da,0, &Mx,0,0, 0,0,0, &dof,0,0,0,0,0);CHKERRQ(ierr); 1123 hx = (ctx->xmax - ctx->xmin)/Mx; 1124 ierr = DMGlobalToLocalBegin(da,X,INSERT_VALUES,Xloc);CHKERRQ(ierr); 1125 ierr = DMGlobalToLocalEnd (da,X,INSERT_VALUES,Xloc);CHKERRQ(ierr); 1126 1127 ierr = VecZeroEntries(F);CHKERRQ(ierr); 1128 1129 ierr = DMDAVecGetArray(da,Xloc,&x);CHKERRQ(ierr); 1130 ierr = DMDAVecGetArray(da,F,&f);CHKERRQ(ierr); 1131 ierr = DMDAGetArray(da,PETSC_TRUE,&slope);CHKERRQ(ierr); 1132 1133 ierr = DMDAGetCorners(da,&xs,0,0,&xm,0,0);CHKERRQ(ierr); 1134 1135 if (ctx->bctype == FVBC_OUTFLOW) { 1136 for (i=xs-2; i<0; i++) { 1137 for (j=0; j<dof; j++) x[i*dof+j] = x[j]; 1138 } 1139 for (i=Mx; i<xs+xm+2; i++) { 1140 for (j=0; j<dof; j++) x[i*dof+j] = x[(xs+xm-1)*dof+j]; 1141 } 1142 } 1143 for (i=xs-1; i<xs+xm+1; i++) { 1144 struct _LimitInfo info; 1145 PetscScalar *cjmpL,*cjmpR; 1146 /* Determine the right eigenvectors R, where A = R \Lambda R^{-1} */ 1147 ierr = (*ctx->physics.characteristic)(ctx->physics.user,dof,&x[i*dof],ctx->R,ctx->Rinv,ctx->speeds);CHKERRQ(ierr); 1148 /* Evaluate jumps across interfaces (i-1, i) and (i, i+1), put in characteristic basis */ 1149 ierr = PetscArrayzero(ctx->cjmpLR,2*dof);CHKERRQ(ierr); 1150 cjmpL = &ctx->cjmpLR[0]; 1151 cjmpR = &ctx->cjmpLR[dof]; 1152 for (j=0; j<dof; j++) { 1153 PetscScalar jmpL,jmpR; 1154 jmpL = x[(i+0)*dof+j] - x[(i-1)*dof+j]; 1155 jmpR = x[(i+1)*dof+j] - x[(i+0)*dof+j]; 1156 for (k=0; k<dof; k++) { 1157 cjmpL[k] += ctx->Rinv[k+j*dof] * jmpL; 1158 cjmpR[k] += ctx->Rinv[k+j*dof] * jmpR; 1159 } 1160 } 1161 /* Apply limiter to the left and right characteristic jumps */ 1162 info.m = dof; 1163 info.hx = hx; 1164 (*ctx->limit)(&info,cjmpL,cjmpR,ctx->cslope); 1165 for (j=0; j<dof; j++) ctx->cslope[j] /= hx; /* rescale to a slope */ 1166 for (j=0; j<dof; j++) { 1167 PetscScalar tmp = 0; 1168 for (k=0; k<dof; k++) tmp += ctx->R[j+k*dof] * ctx->cslope[k]; 1169 slope[i*dof+j] = tmp; 1170 } 1171 } 1172 1173 for (i=xs; i<xs+xm+1; i++) { 1174 PetscReal maxspeed; 1175 PetscScalar *uL,*uR; 1176 uL = &ctx->uLR[0]; 1177 uR = &ctx->uLR[dof]; 1178 for (j=0; j<dof; j++) { 1179 uL[j] = x[(i-1)*dof+j] + slope[(i-1)*dof+j]*hx/2; 1180 uR[j] = x[(i-0)*dof+j] - slope[(i-0)*dof+j]*hx/2; 1181 } 1182 ierr = (*ctx->physics.riemann)(ctx->physics.user,dof,uL,uR,ctx->flux,&maxspeed);CHKERRQ(ierr); 1183 cfl_idt = PetscMax(cfl_idt,PetscAbsScalar(maxspeed/hx)); /* Max allowable value of 1/Delta t */ 1184 1185 if (i > xs) { 1186 for (j=0; j<dof; j++) f[(i-1)*dof+j] -= ctx->flux[j]/hx; 1187 } 1188 if (i < xs+xm) { 1189 for (j=0; j<dof; j++) f[i*dof+j] += ctx->flux[j]/hx; 1190 } 1191 } 1192 1193 ierr = DMDAVecRestoreArray(da,Xloc,&x);CHKERRQ(ierr); 1194 ierr = DMDAVecRestoreArray(da,F,&f);CHKERRQ(ierr); 1195 ierr = DMDARestoreArray(da,PETSC_TRUE,&slope);CHKERRQ(ierr); 1196 ierr = DMRestoreLocalVector(da,&Xloc);CHKERRQ(ierr); 1197 1198 ierr = MPI_Allreduce(&cfl_idt,&ctx->cfl_idt,1,MPIU_REAL,MPIU_MAX,PetscObjectComm((PetscObject)da));CHKERRMPI(ierr); 1199 if (0) { 1200 /* We need to a way to inform the TS of a CFL constraint, this is a debugging fragment */ 1201 PetscReal dt,tnow; 1202 ierr = TSGetTimeStep(ts,&dt);CHKERRQ(ierr); 1203 ierr = TSGetTime(ts,&tnow);CHKERRQ(ierr); 1204 if (dt > 0.5/ctx->cfl_idt) { 1205 ierr = PetscPrintf(ctx->comm,"Stability constraint exceeded at t=%g, dt %g > %g\n",(double)tnow,(double)dt,(double)(0.5/ctx->cfl_idt));CHKERRQ(ierr); 1206 } 1207 } 1208 PetscFunctionReturn(0); 1209 } 1210 1211 static PetscErrorCode SmallMatMultADB(PetscScalar *C,PetscInt bs,const PetscScalar *A,const PetscReal *D,const PetscScalar *B) 1212 { 1213 PetscInt i,j,k; 1214 1215 PetscFunctionBeginUser; 1216 for (i=0; i<bs; i++) { 1217 for (j=0; j<bs; j++) { 1218 PetscScalar tmp = 0; 1219 for (k=0; k<bs; k++) tmp += A[i*bs+k] * D[k] * B[k*bs+j]; 1220 C[i*bs+j] = tmp; 1221 } 1222 } 1223 PetscFunctionReturn(0); 1224 } 1225 1226 static PetscErrorCode FVIJacobian(TS ts,PetscReal t,Vec X,Vec Xdot,PetscReal shift,Mat A,Mat B,void *vctx) 1227 { 1228 FVCtx *ctx = (FVCtx*)vctx; 1229 PetscErrorCode ierr; 1230 PetscInt i,j,dof = ctx->physics.dof; 1231 PetscScalar *J; 1232 const PetscScalar *x; 1233 PetscReal hx; 1234 DM da; 1235 DMDALocalInfo dainfo; 1236 1237 PetscFunctionBeginUser; 1238 ierr = TSGetDM(ts,&da);CHKERRQ(ierr); 1239 ierr = DMDAVecGetArrayRead(da,X,(void*)&x);CHKERRQ(ierr); 1240 ierr = DMDAGetLocalInfo(da,&dainfo);CHKERRQ(ierr); 1241 hx = (ctx->xmax - ctx->xmin)/dainfo.mx; 1242 ierr = PetscMalloc1(dof*dof,&J);CHKERRQ(ierr); 1243 for (i=dainfo.xs; i<dainfo.xs+dainfo.xm; i++) { 1244 ierr = (*ctx->physics.characteristic)(ctx->physics.user,dof,&x[i*dof],ctx->R,ctx->Rinv,ctx->speeds);CHKERRQ(ierr); 1245 for (j=0; j<dof; j++) ctx->speeds[j] = PetscAbs(ctx->speeds[j]); 1246 ierr = SmallMatMultADB(J,dof,ctx->R,ctx->speeds,ctx->Rinv);CHKERRQ(ierr); 1247 for (j=0; j<dof*dof; j++) J[j] = J[j]/hx + shift*(j/dof == j%dof); 1248 ierr = MatSetValuesBlocked(B,1,&i,1,&i,J,INSERT_VALUES);CHKERRQ(ierr); 1249 } 1250 ierr = PetscFree(J);CHKERRQ(ierr); 1251 ierr = DMDAVecRestoreArrayRead(da,X,(void*)&x);CHKERRQ(ierr); 1252 1253 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1254 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1255 if (A != B) { 1256 ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1257 ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1258 } 1259 PetscFunctionReturn(0); 1260 } 1261 1262 static PetscErrorCode FVSample(FVCtx *ctx,DM da,PetscReal time,Vec U) 1263 { 1264 PetscErrorCode ierr; 1265 PetscScalar *u,*uj; 1266 PetscInt i,j,k,dof,xs,xm,Mx; 1267 1268 PetscFunctionBeginUser; 1269 if (!ctx->physics.sample) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Physics has not provided a sampling function"); 1270 ierr = DMDAGetInfo(da,0, &Mx,0,0, 0,0,0, &dof,0,0,0,0,0);CHKERRQ(ierr); 1271 ierr = DMDAGetCorners(da,&xs,0,0,&xm,0,0);CHKERRQ(ierr); 1272 ierr = DMDAVecGetArray(da,U,&u);CHKERRQ(ierr); 1273 ierr = PetscMalloc1(dof,&uj);CHKERRQ(ierr); 1274 for (i=xs; i<xs+xm; i++) { 1275 const PetscReal h = (ctx->xmax-ctx->xmin)/Mx,xi = ctx->xmin+h/2+i*h; 1276 const PetscInt N = 200; 1277 /* Integrate over cell i using trapezoid rule with N points. */ 1278 for (k=0; k<dof; k++) u[i*dof+k] = 0; 1279 for (j=0; j<N+1; j++) { 1280 PetscScalar xj = xi+h*(j-N/2)/(PetscReal)N; 1281 ierr = (*ctx->physics.sample)(ctx->physics.user,ctx->initial,ctx->bctype,ctx->xmin,ctx->xmax,time,xj,uj);CHKERRQ(ierr); 1282 for (k=0; k<dof; k++) u[i*dof+k] += ((j==0 || j==N) ? 0.5 : 1.0)*uj[k]/N; 1283 } 1284 } 1285 ierr = DMDAVecRestoreArray(da,U,&u);CHKERRQ(ierr); 1286 ierr = PetscFree(uj);CHKERRQ(ierr); 1287 PetscFunctionReturn(0); 1288 } 1289 1290 static PetscErrorCode SolutionStatsView(DM da,Vec X,PetscViewer viewer) 1291 { 1292 PetscErrorCode ierr; 1293 PetscReal xmin,xmax; 1294 PetscScalar sum,tvsum,tvgsum; 1295 const PetscScalar *x; 1296 PetscInt imin,imax,Mx,i,j,xs,xm,dof; 1297 Vec Xloc; 1298 PetscBool iascii; 1299 1300 PetscFunctionBeginUser; 1301 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 1302 if (iascii) { 1303 /* PETSc lacks a function to compute total variation norm (difficult in multiple dimensions), we do it here */ 1304 ierr = DMGetLocalVector(da,&Xloc);CHKERRQ(ierr); 1305 ierr = DMGlobalToLocalBegin(da,X,INSERT_VALUES,Xloc);CHKERRQ(ierr); 1306 ierr = DMGlobalToLocalEnd (da,X,INSERT_VALUES,Xloc);CHKERRQ(ierr); 1307 ierr = DMDAVecGetArrayRead(da,Xloc,(void*)&x);CHKERRQ(ierr); 1308 ierr = DMDAGetCorners(da,&xs,0,0,&xm,0,0);CHKERRQ(ierr); 1309 ierr = DMDAGetInfo(da,0, &Mx,0,0, 0,0,0, &dof,0,0,0,0,0);CHKERRQ(ierr); 1310 tvsum = 0; 1311 for (i=xs; i<xs+xm; i++) { 1312 for (j=0; j<dof; j++) tvsum += PetscAbsScalar(x[i*dof+j] - x[(i-1)*dof+j]); 1313 } 1314 ierr = MPI_Allreduce(&tvsum,&tvgsum,1,MPIU_REAL,MPIU_SUM,PetscObjectComm((PetscObject)da));CHKERRMPI(ierr); 1315 ierr = DMDAVecRestoreArrayRead(da,Xloc,(void*)&x);CHKERRQ(ierr); 1316 ierr = DMRestoreLocalVector(da,&Xloc);CHKERRQ(ierr); 1317 1318 ierr = VecMin(X,&imin,&xmin);CHKERRQ(ierr); 1319 ierr = VecMax(X,&imax,&xmax);CHKERRQ(ierr); 1320 ierr = VecSum(X,&sum);CHKERRQ(ierr); 1321 ierr = PetscViewerASCIIPrintf(viewer,"Solution range [%8.5f,%8.5f] with extrema at %D and %D, mean %8.5f, ||x||_TV %8.5f\n",(double)xmin,(double)xmax,imin,imax,(double)(sum/Mx),(double)(tvgsum/Mx));CHKERRQ(ierr); 1322 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Viewer type not supported"); 1323 PetscFunctionReturn(0); 1324 } 1325 1326 static PetscErrorCode SolutionErrorNorms(FVCtx *ctx,DM da,PetscReal t,Vec X,PetscReal *nrm1,PetscReal *nrmsup) 1327 { 1328 PetscErrorCode ierr; 1329 Vec Y; 1330 PetscInt Mx; 1331 1332 PetscFunctionBeginUser; 1333 ierr = VecGetSize(X,&Mx);CHKERRQ(ierr); 1334 ierr = VecDuplicate(X,&Y);CHKERRQ(ierr); 1335 ierr = FVSample(ctx,da,t,Y);CHKERRQ(ierr); 1336 ierr = VecAYPX(Y,-1,X);CHKERRQ(ierr); 1337 ierr = VecNorm(Y,NORM_1,nrm1);CHKERRQ(ierr); 1338 ierr = VecNorm(Y,NORM_INFINITY,nrmsup);CHKERRQ(ierr); 1339 *nrm1 /= Mx; 1340 ierr = VecDestroy(&Y);CHKERRQ(ierr); 1341 PetscFunctionReturn(0); 1342 } 1343 1344 int main(int argc,char *argv[]) 1345 { 1346 char lname[256] = "mc",physname[256] = "advect",final_fname[256] = "solution.m"; 1347 PetscFunctionList limiters = 0,physics = 0; 1348 MPI_Comm comm; 1349 TS ts; 1350 DM da; 1351 Vec X,X0,R; 1352 Mat B; 1353 FVCtx ctx; 1354 PetscInt i,dof,xs,xm,Mx,draw = 0; 1355 PetscBool view_final = PETSC_FALSE; 1356 PetscReal ptime; 1357 PetscErrorCode ierr; 1358 1359 ierr = PetscInitialize(&argc,&argv,0,help);if (ierr) return ierr; 1360 comm = PETSC_COMM_WORLD; 1361 ierr = PetscMemzero(&ctx,sizeof(ctx));CHKERRQ(ierr); 1362 1363 /* Register limiters to be available on the command line */ 1364 ierr = PetscFunctionListAdd(&limiters,"upwind" ,Limit_Upwind);CHKERRQ(ierr); 1365 ierr = PetscFunctionListAdd(&limiters,"lax-wendroff" ,Limit_LaxWendroff);CHKERRQ(ierr); 1366 ierr = PetscFunctionListAdd(&limiters,"beam-warming" ,Limit_BeamWarming);CHKERRQ(ierr); 1367 ierr = PetscFunctionListAdd(&limiters,"fromm" ,Limit_Fromm);CHKERRQ(ierr); 1368 ierr = PetscFunctionListAdd(&limiters,"minmod" ,Limit_Minmod);CHKERRQ(ierr); 1369 ierr = PetscFunctionListAdd(&limiters,"superbee" ,Limit_Superbee);CHKERRQ(ierr); 1370 ierr = PetscFunctionListAdd(&limiters,"mc" ,Limit_MC);CHKERRQ(ierr); 1371 ierr = PetscFunctionListAdd(&limiters,"vanleer" ,Limit_VanLeer);CHKERRQ(ierr); 1372 ierr = PetscFunctionListAdd(&limiters,"vanalbada" ,Limit_VanAlbada);CHKERRQ(ierr); 1373 ierr = PetscFunctionListAdd(&limiters,"vanalbadatvd" ,Limit_VanAlbadaTVD);CHKERRQ(ierr); 1374 ierr = PetscFunctionListAdd(&limiters,"koren" ,Limit_Koren);CHKERRQ(ierr); 1375 ierr = PetscFunctionListAdd(&limiters,"korensym" ,Limit_KorenSym);CHKERRQ(ierr); 1376 ierr = PetscFunctionListAdd(&limiters,"koren3" ,Limit_Koren3);CHKERRQ(ierr); 1377 ierr = PetscFunctionListAdd(&limiters,"cada-torrilhon2" ,Limit_CadaTorrilhon2);CHKERRQ(ierr); 1378 ierr = PetscFunctionListAdd(&limiters,"cada-torrilhon3-r0p1",Limit_CadaTorrilhon3R0p1);CHKERRQ(ierr); 1379 ierr = PetscFunctionListAdd(&limiters,"cada-torrilhon3-r1" ,Limit_CadaTorrilhon3R1);CHKERRQ(ierr); 1380 ierr = PetscFunctionListAdd(&limiters,"cada-torrilhon3-r10" ,Limit_CadaTorrilhon3R10);CHKERRQ(ierr); 1381 ierr = PetscFunctionListAdd(&limiters,"cada-torrilhon3-r100",Limit_CadaTorrilhon3R100);CHKERRQ(ierr); 1382 1383 /* Register physical models to be available on the command line */ 1384 ierr = PetscFunctionListAdd(&physics,"advect" ,PhysicsCreate_Advect);CHKERRQ(ierr); 1385 ierr = PetscFunctionListAdd(&physics,"burgers" ,PhysicsCreate_Burgers);CHKERRQ(ierr); 1386 ierr = PetscFunctionListAdd(&physics,"traffic" ,PhysicsCreate_Traffic);CHKERRQ(ierr); 1387 ierr = PetscFunctionListAdd(&physics,"acoustics" ,PhysicsCreate_Acoustics);CHKERRQ(ierr); 1388 ierr = PetscFunctionListAdd(&physics,"isogas" ,PhysicsCreate_IsoGas);CHKERRQ(ierr); 1389 ierr = PetscFunctionListAdd(&physics,"shallow" ,PhysicsCreate_Shallow);CHKERRQ(ierr); 1390 1391 ctx.comm = comm; 1392 ctx.cfl = 0.9; ctx.bctype = FVBC_PERIODIC; 1393 ctx.xmin = -1; ctx.xmax = 1; 1394 ierr = PetscOptionsBegin(comm,NULL,"Finite Volume solver options","");CHKERRQ(ierr); 1395 ierr = PetscOptionsReal("-xmin","X min","",ctx.xmin,&ctx.xmin,NULL);CHKERRQ(ierr); 1396 ierr = PetscOptionsReal("-xmax","X max","",ctx.xmax,&ctx.xmax,NULL);CHKERRQ(ierr); 1397 ierr = PetscOptionsFList("-limit","Name of flux limiter to use","",limiters,lname,lname,sizeof(lname),NULL);CHKERRQ(ierr); 1398 ierr = PetscOptionsFList("-physics","Name of physics (Riemann solver and characteristics) to use","",physics,physname,physname,sizeof(physname),NULL);CHKERRQ(ierr); 1399 ierr = PetscOptionsInt("-draw","Draw solution vector, bitwise OR of (1=initial,2=final,4=final error)","",draw,&draw,NULL);CHKERRQ(ierr); 1400 ierr = PetscOptionsString("-view_final","Write final solution in ASCII MATLAB format to given file name","",final_fname,final_fname,sizeof(final_fname),&view_final);CHKERRQ(ierr); 1401 ierr = PetscOptionsInt("-initial","Initial condition (depends on the physics)","",ctx.initial,&ctx.initial,NULL);CHKERRQ(ierr); 1402 ierr = PetscOptionsBool("-exact","Compare errors with exact solution","",ctx.exact,&ctx.exact,NULL);CHKERRQ(ierr); 1403 ierr = PetscOptionsReal("-cfl","CFL number to time step at","",ctx.cfl,&ctx.cfl,NULL);CHKERRQ(ierr); 1404 ierr = PetscOptionsEnum("-bc_type","Boundary condition","",FVBCTypes,(PetscEnum)ctx.bctype,(PetscEnum*)&ctx.bctype,NULL);CHKERRQ(ierr); 1405 ierr = PetscOptionsEnd();CHKERRQ(ierr); 1406 1407 /* Choose the limiter from the list of registered limiters */ 1408 ierr = PetscFunctionListFind(limiters,lname,&ctx.limit);CHKERRQ(ierr); 1409 if (!ctx.limit) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Limiter '%s' not found",lname); 1410 1411 /* Choose the physics from the list of registered models */ 1412 { 1413 PetscErrorCode (*r)(FVCtx*); 1414 ierr = PetscFunctionListFind(physics,physname,&r);CHKERRQ(ierr); 1415 if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Physics '%s' not found",physname); 1416 /* Create the physics, will set the number of fields and their names */ 1417 ierr = (*r)(&ctx);CHKERRQ(ierr); 1418 } 1419 1420 /* Create a DMDA to manage the parallel grid */ 1421 ierr = DMDACreate1d(comm,DM_BOUNDARY_PERIODIC,50,ctx.physics.dof,2,NULL,&da);CHKERRQ(ierr); 1422 ierr = DMSetFromOptions(da);CHKERRQ(ierr); 1423 ierr = DMSetUp(da);CHKERRQ(ierr); 1424 /* Inform the DMDA of the field names provided by the physics. */ 1425 /* The names will be shown in the title bars when run with -ts_monitor_draw_solution */ 1426 for (i=0; i<ctx.physics.dof; i++) { 1427 ierr = DMDASetFieldName(da,i,ctx.physics.fieldname[i]);CHKERRQ(ierr); 1428 } 1429 ierr = DMDAGetInfo(da,0, &Mx,0,0, 0,0,0, &dof,0,0,0,0,0);CHKERRQ(ierr); 1430 ierr = DMDAGetCorners(da,&xs,0,0,&xm,0,0);CHKERRQ(ierr); 1431 1432 /* Set coordinates of cell centers */ 1433 ierr = DMDASetUniformCoordinates(da,ctx.xmin+0.5*(ctx.xmax-ctx.xmin)/Mx,ctx.xmax+0.5*(ctx.xmax-ctx.xmin)/Mx,0,0,0,0);CHKERRQ(ierr); 1434 1435 /* Allocate work space for the Finite Volume solver (so it doesn't have to be reallocated on each function evaluation) */ 1436 ierr = PetscMalloc4(dof*dof,&ctx.R,dof*dof,&ctx.Rinv,2*dof,&ctx.cjmpLR,1*dof,&ctx.cslope);CHKERRQ(ierr); 1437 ierr = PetscMalloc3(2*dof,&ctx.uLR,dof,&ctx.flux,dof,&ctx.speeds);CHKERRQ(ierr); 1438 1439 /* Create a vector to store the solution and to save the initial state */ 1440 ierr = DMCreateGlobalVector(da,&X);CHKERRQ(ierr); 1441 ierr = VecDuplicate(X,&X0);CHKERRQ(ierr); 1442 ierr = VecDuplicate(X,&R);CHKERRQ(ierr); 1443 1444 ierr = DMCreateMatrix(da,&B);CHKERRQ(ierr); 1445 1446 /* Create a time-stepping object */ 1447 ierr = TSCreate(comm,&ts);CHKERRQ(ierr); 1448 ierr = TSSetDM(ts,da);CHKERRQ(ierr); 1449 ierr = TSSetRHSFunction(ts,R,FVRHSFunction,&ctx);CHKERRQ(ierr); 1450 ierr = TSSetIJacobian(ts,B,B,FVIJacobian,&ctx);CHKERRQ(ierr); 1451 ierr = TSSetType(ts,TSSSP);CHKERRQ(ierr); 1452 ierr = TSSetMaxTime(ts,10);CHKERRQ(ierr); 1453 ierr = TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);CHKERRQ(ierr); 1454 1455 /* Compute initial conditions and starting time step */ 1456 ierr = FVSample(&ctx,da,0,X0);CHKERRQ(ierr); 1457 ierr = FVRHSFunction(ts,0,X0,X,(void*)&ctx);CHKERRQ(ierr); /* Initial function evaluation, only used to determine max speed */ 1458 ierr = VecCopy(X0,X);CHKERRQ(ierr); /* The function value was not used so we set X=X0 again */ 1459 ierr = TSSetTimeStep(ts,ctx.cfl/ctx.cfl_idt);CHKERRQ(ierr); 1460 ierr = TSSetFromOptions(ts);CHKERRQ(ierr); /* Take runtime options */ 1461 ierr = SolutionStatsView(da,X,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 1462 { 1463 PetscReal nrm1,nrmsup; 1464 PetscInt steps; 1465 1466 ierr = TSSolve(ts,X);CHKERRQ(ierr); 1467 ierr = TSGetSolveTime(ts,&ptime);CHKERRQ(ierr); 1468 ierr = TSGetStepNumber(ts,&steps);CHKERRQ(ierr); 1469 1470 ierr = PetscPrintf(comm,"Final time %8.5f, steps %D\n",(double)ptime,steps);CHKERRQ(ierr); 1471 if (ctx.exact) { 1472 ierr = SolutionErrorNorms(&ctx,da,ptime,X,&nrm1,&nrmsup);CHKERRQ(ierr); 1473 ierr = PetscPrintf(comm,"Error ||x-x_e||_1 %8.4e ||x-x_e||_sup %8.4e\n",(double)nrm1,(double)nrmsup);CHKERRQ(ierr); 1474 } 1475 } 1476 1477 ierr = SolutionStatsView(da,X,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 1478 if (draw & 0x1) {ierr = VecView(X0,PETSC_VIEWER_DRAW_WORLD);CHKERRQ(ierr);} 1479 if (draw & 0x2) {ierr = VecView(X,PETSC_VIEWER_DRAW_WORLD);CHKERRQ(ierr);} 1480 if (draw & 0x4) { 1481 Vec Y; 1482 ierr = VecDuplicate(X,&Y);CHKERRQ(ierr); 1483 ierr = FVSample(&ctx,da,ptime,Y);CHKERRQ(ierr); 1484 ierr = VecAYPX(Y,-1,X);CHKERRQ(ierr); 1485 ierr = VecView(Y,PETSC_VIEWER_DRAW_WORLD);CHKERRQ(ierr); 1486 ierr = VecDestroy(&Y);CHKERRQ(ierr); 1487 } 1488 1489 if (view_final) { 1490 PetscViewer viewer; 1491 ierr = PetscViewerASCIIOpen(PETSC_COMM_WORLD,final_fname,&viewer);CHKERRQ(ierr); 1492 ierr = PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_MATLAB);CHKERRQ(ierr); 1493 ierr = VecView(X,viewer);CHKERRQ(ierr); 1494 ierr = PetscViewerPopFormat(viewer);CHKERRQ(ierr); 1495 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 1496 } 1497 1498 /* Clean up */ 1499 ierr = (*ctx.physics.destroy)(ctx.physics.user);CHKERRQ(ierr); 1500 for (i=0; i<ctx.physics.dof; i++) {ierr = PetscFree(ctx.physics.fieldname[i]);CHKERRQ(ierr);} 1501 ierr = PetscFree4(ctx.R,ctx.Rinv,ctx.cjmpLR,ctx.cslope);CHKERRQ(ierr); 1502 ierr = PetscFree3(ctx.uLR,ctx.flux,ctx.speeds);CHKERRQ(ierr); 1503 ierr = VecDestroy(&X);CHKERRQ(ierr); 1504 ierr = VecDestroy(&X0);CHKERRQ(ierr); 1505 ierr = VecDestroy(&R);CHKERRQ(ierr); 1506 ierr = MatDestroy(&B);CHKERRQ(ierr); 1507 ierr = DMDestroy(&da);CHKERRQ(ierr); 1508 ierr = TSDestroy(&ts);CHKERRQ(ierr); 1509 ierr = PetscFunctionListDestroy(&limiters);CHKERRQ(ierr); 1510 ierr = PetscFunctionListDestroy(&physics);CHKERRQ(ierr); 1511 ierr = PetscFinalize(); 1512 return ierr; 1513 } 1514 1515 /*TEST 1516 1517 build: 1518 requires: !complex 1519 1520 test: 1521 args: -da_grid_x 100 -initial 1 -xmin -2 -xmax 5 -exact -limit mc 1522 requires: !complex !single 1523 1524 test: 1525 suffix: 2 1526 args: -da_grid_x 100 -initial 2 -xmin -2 -xmax 2 -exact -limit mc -physics burgers -bc_type outflow -ts_max_time 1 1527 filter: sed "s/at 48/at 0/g" 1528 requires: !complex !single 1529 1530 test: 1531 suffix: 3 1532 args: -da_grid_x 100 -initial 2 -xmin -2 -xmax 2 -exact -limit mc -physics burgers -bc_type outflow -ts_max_time 1 1533 nsize: 3 1534 filter: sed "s/at 48/at 0/g" 1535 requires: !complex !single 1536 1537 TEST*/ 1538