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