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