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