1 /* 2 Note: 3 -hratio is the ratio between mesh size of coarse grids and fine grids 4 */ 5 6 static const char help[] = "1D periodic Finite Volume solver in slope-limiter form with semidiscrete time stepping.\n" 7 " advect - Constant coefficient scalar advection\n" 8 " u_t + (a*u)_x = 0\n" 9 " shallow - 1D Shallow water equations (Saint Venant System)\n" 10 " h_t + (q)_x = 0 \n" 11 " q_t + (\frac{q^2}{h} + g/2*h^2)_x = 0 \n" 12 " where, h(x,t) denotes the height of the water and q(x,t) the momentum.\n" 13 " for this toy problem, we choose different meshsizes for different sub-domains (slow-fast-slow), say\n" 14 " hxs = hratio*hxf \n" 15 " where hxs and hxf are the grid spacings for coarse and fine grids respectively.\n" 16 " exact - Exact Riemann solver which usually needs to perform a Newton iteration to connect\n" 17 " the states across shocks and rarefactions\n" 18 " simulation - use reference solution which is generated by smaller time step size to be true solution,\n" 19 " also the reference solution should be generated by user and stored in a binary file.\n" 20 " characteristic - Limit the characteristic variables, this is usually preferred (default)\n" 21 " bc_type - Boundary condition for the problem, options are: periodic, outflow, inflow " 22 "Several problem descriptions (initial data, physics specific features, boundary data) can be chosen with -initial N\n\n" 23 "The problem size should be set with -da_grid_x M\n\n"; 24 25 /* 26 Example: 27 ./ex4 -da_grid_x 40 -initial 1 -hratio 1 -limit mc -ts_dt 0.01 -ts_max_time 7.0 -ts_type mprk -ts_mprk_type 2a22 -ts_monitor_draw_solution -physics shallow -bc_type outflow -xmin 0 -xmax 50 -ts_use_splitrhsfunction 0 28 ./ex4 -da_grid_x 40 -initial 2 -hratio 1 -limit mc -ts_dt 0.01 -ts_max_time 2.5 -ts_type mprk -ts_mprk_type 2a22 -ts_monitor_draw_solution -physics shallow -bc_type outflow -xmin 0 -xmax 50 -ts_use_splitrhsfunction 0 29 ./ex4 -da_grid_x 40 -initial 3 -hratio 1 -limit mc -ts_dt 0.01 -ts_max_time 4.0 -ts_type mprk -ts_mprk_type 2a22 -ts_monitor_draw_solution -physics shallow -bc_type outflow -xmin 0 -xmax 50 -ts_use_splitrhsfunction 0 30 ./ex4 -da_grid_x 40 -initial 4 -hratio 1 -limit koren3 -ts_dt 0.01 -ts_max_time 4.0 -ts_type mprk -ts_mprk_type 2a22 -ts_monitor_draw_solution -physics shallow -bc_type outflow -xmin 0 -xmax 50 -ts_use_splitrhsfunction 0 31 ./ex4 -da_grid_x 40 -initial 5 -hratio 1 -limit mc -ts_dt 0.01 -ts_max_time 5.0 -ts_type mprk -ts_mprk_type 2a22 -ts_monitor_draw_solution -physics shallow -bc_type outflow -xmin 0 -xmax 50 -ts_use_splitrhsfunction 0 32 33 Contributed by: Aidan Hamilton <aidan@udel.edu> 34 */ 35 36 #include <petscts.h> 37 #include <petscdm.h> 38 #include <petscdmda.h> 39 #include <petscdraw.h> 40 #include "finitevolume1d.h" 41 #include <petsc/private/kernels/blockinvert.h> 42 43 static inline PetscReal RangeMod(PetscReal a, PetscReal xmin, PetscReal xmax) 44 { 45 PetscReal range = xmax - xmin; 46 return xmin + PetscFmodReal(range + PetscFmodReal(a, range), range); 47 } 48 static inline PetscReal MaxAbs(PetscReal a, PetscReal b) 49 { 50 return (PetscAbs(a) > PetscAbs(b)) ? a : b; 51 } 52 /* --------------------------------- Advection ----------------------------------- */ 53 typedef struct { 54 PetscReal a; /* advective velocity */ 55 } AdvectCtx; 56 57 static PetscErrorCode PhysicsRiemann_Advect(void *vctx, PetscInt m, const PetscScalar *uL, const PetscScalar *uR, PetscScalar *flux, PetscReal *maxspeed) 58 { 59 AdvectCtx *ctx = (AdvectCtx *)vctx; 60 PetscReal speed; 61 62 PetscFunctionBeginUser; 63 speed = ctx->a; 64 flux[0] = PetscMax(0, speed) * uL[0] + PetscMin(0, speed) * uR[0]; 65 *maxspeed = speed; 66 PetscFunctionReturn(PETSC_SUCCESS); 67 } 68 69 static PetscErrorCode PhysicsCharacteristic_Advect(void *vctx, PetscInt m, const PetscScalar *u, PetscScalar *X, PetscScalar *Xi, PetscReal *speeds) 70 { 71 AdvectCtx *ctx = (AdvectCtx *)vctx; 72 73 PetscFunctionBeginUser; 74 X[0] = 1.; 75 Xi[0] = 1.; 76 speeds[0] = ctx->a; 77 PetscFunctionReturn(PETSC_SUCCESS); 78 } 79 80 static PetscErrorCode PhysicsSample_Advect(void *vctx, PetscInt initial, FVBCType bctype, PetscReal xmin, PetscReal xmax, PetscReal t, PetscReal x, PetscReal *u) 81 { 82 AdvectCtx *ctx = (AdvectCtx *)vctx; 83 PetscReal a = ctx->a, x0; 84 85 PetscFunctionBeginUser; 86 switch (bctype) { 87 case FVBC_OUTFLOW: 88 x0 = x - a * t; 89 break; 90 case FVBC_PERIODIC: 91 x0 = RangeMod(x - a * t, xmin, xmax); 92 break; 93 default: 94 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "unknown BCType"); 95 } 96 switch (initial) { 97 case 0: 98 u[0] = (x0 < 0) ? 1 : -1; 99 break; 100 case 1: 101 u[0] = (x0 < 0) ? -1 : 1; 102 break; 103 case 2: 104 u[0] = (0 < x0 && x0 < 1) ? 1 : 0; 105 break; 106 case 3: 107 u[0] = PetscSinReal(2 * PETSC_PI * x0); 108 break; 109 case 4: 110 u[0] = PetscAbs(x0); 111 break; 112 case 5: 113 u[0] = (x0 < 0 || x0 > 0.5) ? 0 : PetscSqr(PetscSinReal(2 * PETSC_PI * x0)); 114 break; 115 case 6: 116 u[0] = (x0 < 0) ? 0 : ((x0 < 1) ? x0 : ((x0 < 2) ? 2 - x0 : 0)); 117 break; 118 case 7: 119 u[0] = PetscPowReal(PetscSinReal(PETSC_PI * x0), 10.0); 120 break; 121 default: 122 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "unknown initial condition"); 123 } 124 PetscFunctionReturn(PETSC_SUCCESS); 125 } 126 127 static PetscErrorCode PhysicsCreate_Advect(FVCtx *ctx) 128 { 129 AdvectCtx *user; 130 131 PetscFunctionBeginUser; 132 PetscCall(PetscNew(&user)); 133 ctx->physics2.sample2 = PhysicsSample_Advect; 134 ctx->physics2.riemann2 = PhysicsRiemann_Advect; 135 ctx->physics2.characteristic2 = PhysicsCharacteristic_Advect; 136 ctx->physics2.destroy = PhysicsDestroy_SimpleFree; 137 ctx->physics2.user = user; 138 ctx->physics2.dof = 1; 139 140 PetscCall(PetscStrallocpy("u", &ctx->physics2.fieldname[0])); 141 user->a = 1; 142 PetscOptionsBegin(ctx->comm, ctx->prefix, "Options for advection", ""); 143 PetscCall(PetscOptionsReal("-physics_advect_a", "Speed", "", user->a, &user->a, NULL)); 144 PetscOptionsEnd(); 145 PetscFunctionReturn(PETSC_SUCCESS); 146 } 147 /* --------------------------------- Shallow Water ----------------------------------- */ 148 149 typedef struct { 150 PetscReal gravity; 151 } ShallowCtx; 152 153 static inline void ShallowFlux(ShallowCtx *phys, const PetscScalar *u, PetscScalar *f) 154 { 155 f[0] = u[1]; 156 f[1] = PetscSqr(u[1]) / u[0] + 0.5 * phys->gravity * PetscSqr(u[0]); 157 } 158 159 static PetscErrorCode PhysicsRiemann_Shallow_Exact(void *vctx, PetscInt m, const PetscScalar *uL, const PetscScalar *uR, PetscScalar *flux, PetscReal *maxspeed) 160 { 161 ShallowCtx *phys = (ShallowCtx *)vctx; 162 PetscScalar g = phys->gravity, ustar[2], cL, cR, c, cstar; 163 struct { 164 PetscScalar h, u; 165 } L = {uL[0], uL[1] / uL[0]}, R = {uR[0], uR[1] / uR[0]}, star; 166 PetscInt i; 167 168 PetscFunctionBeginUser; 169 PetscCheck((L.h > 0 && R.h > 0), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Reconstructed thickness is negative"); 170 cL = PetscSqrtScalar(g * L.h); 171 cR = PetscSqrtScalar(g * R.h); 172 c = PetscMax(cL, cR); 173 { 174 /* Solve for star state */ 175 const PetscInt maxits = 50; 176 PetscScalar tmp, res, res0 = 0, h0, h = 0.5 * (L.h + R.h); /* initial guess */ 177 h0 = h; 178 for (i = 0; i < maxits; i++) { 179 PetscScalar fr, fl, dfr, dfl; 180 fl = (L.h < h) ? PetscSqrtScalar(0.5 * g * (h * h - L.h * L.h) * (1 / L.h - 1 / h)) /* shock */ 181 : 2 * PetscSqrtScalar(g * h) - 2 * PetscSqrtScalar(g * L.h); /* rarefaction */ 182 fr = (R.h < h) ? PetscSqrtScalar(0.5 * g * (h * h - R.h * R.h) * (1 / R.h - 1 / h)) /* shock */ 183 : 2 * PetscSqrtScalar(g * h) - 2 * PetscSqrtScalar(g * R.h); /* rarefaction */ 184 res = R.u - L.u + fr + fl; 185 PetscCheck(!PetscIsInfOrNanScalar(res), PETSC_COMM_SELF, PETSC_ERR_FP, "Infinity or Not-a-Number generated in computation"); 186 if (PetscAbsScalar(res) < PETSC_SQRT_MACHINE_EPSILON || (i > 0 && PetscAbsScalar(h - h0) < PETSC_SQRT_MACHINE_EPSILON)) { 187 star.h = h; 188 star.u = L.u - fl; 189 goto converged; 190 } else if (i > 0 && PetscAbsScalar(res) >= PetscAbsScalar(res0)) { /* Line search */ 191 h = 0.8 * h0 + 0.2 * h; 192 continue; 193 } 194 /* Accept the last step and take another */ 195 res0 = res; 196 h0 = h; 197 dfl = (L.h < h) ? 0.5 / fl * 0.5 * g * (-L.h * L.h / (h * h) - 1 + 2 * h / L.h) : PetscSqrtScalar(g / h); 198 dfr = (R.h < h) ? 0.5 / fr * 0.5 * g * (-R.h * R.h / (h * h) - 1 + 2 * h / R.h) : PetscSqrtScalar(g / h); 199 tmp = h - res / (dfr + dfl); 200 if (tmp <= 0) h /= 2; /* Guard against Newton shooting off to a negative thickness */ 201 else h = tmp; 202 PetscCheck(((h > 0) && PetscIsNormalScalar(h)), PETSC_COMM_SELF, PETSC_ERR_FP, "non-normal iterate h=%g", (double)h); 203 } 204 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_NOT_CONVERGED, "Newton iteration for star.h diverged after %" PetscInt_FMT " iterations", i); 205 } 206 converged: 207 cstar = PetscSqrtScalar(g * star.h); 208 if (L.u - cL < 0 && 0 < star.u - cstar) { /* 1-wave is sonic rarefaction */ 209 PetscScalar ufan[2]; 210 ufan[0] = 1 / g * PetscSqr(L.u / 3 + 2. / 3 * cL); 211 ufan[1] = PetscSqrtScalar(g * ufan[0]) * ufan[0]; 212 ShallowFlux(phys, ufan, flux); 213 } else if (star.u + cstar < 0 && 0 < R.u + cR) { /* 2-wave is sonic rarefaction */ 214 PetscScalar ufan[2]; 215 ufan[0] = 1 / g * PetscSqr(R.u / 3 - 2. / 3 * cR); 216 ufan[1] = -PetscSqrtScalar(g * ufan[0]) * ufan[0]; 217 ShallowFlux(phys, ufan, flux); 218 } 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)) { 219 /* 1-wave is right-travelling shock (supersonic) */ 220 ShallowFlux(phys, uL, flux); 221 } 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)) { 222 /* 2-wave is left-travelling shock (supersonic) */ 223 ShallowFlux(phys, uR, flux); 224 } else { 225 ustar[0] = star.h; 226 ustar[1] = star.h * star.u; 227 ShallowFlux(phys, ustar, flux); 228 } 229 *maxspeed = MaxAbs(MaxAbs(star.u - cstar, star.u + cstar), MaxAbs(L.u - cL, R.u + cR)); 230 PetscFunctionReturn(PETSC_SUCCESS); 231 } 232 233 static PetscErrorCode PhysicsRiemann_Shallow_Rusanov(void *vctx, PetscInt m, const PetscScalar *uL, const PetscScalar *uR, PetscScalar *flux, PetscReal *maxspeed) 234 { 235 ShallowCtx *phys = (ShallowCtx *)vctx; 236 PetscScalar g = phys->gravity, fL[2], fR[2], s; 237 struct { 238 PetscScalar h, u; 239 } L = {uL[0], uL[1] / uL[0]}, R = {uR[0], uR[1] / uR[0]}; 240 PetscReal tol = 1e-6; 241 242 PetscFunctionBeginUser; 243 /* Positivity preserving modification*/ 244 if (L.h < tol) L.u = 0.0; 245 if (R.h < tol) R.u = 0.0; 246 247 /*simple pos preserve limiter*/ 248 if (L.h < 0) L.h = 0; 249 if (R.h < 0) R.h = 0; 250 251 ShallowFlux(phys, uL, fL); 252 ShallowFlux(phys, uR, fR); 253 254 s = PetscMax(PetscAbs(L.u) + PetscSqrtScalar(g * L.h), PetscAbs(R.u) + PetscSqrtScalar(g * R.h)); 255 flux[0] = 0.5 * (fL[0] + fR[0]) + 0.5 * s * (L.h - R.h); 256 flux[1] = 0.5 * (fL[1] + fR[1]) + 0.5 * s * (uL[1] - uR[1]); 257 *maxspeed = s; 258 PetscFunctionReturn(PETSC_SUCCESS); 259 } 260 261 static PetscErrorCode PhysicsCharacteristic_Conservative(void *vctx, PetscInt m, const PetscScalar *u, PetscScalar *X, PetscScalar *Xi, PetscReal *speeds) 262 { 263 PetscInt i, j; 264 265 PetscFunctionBeginUser; 266 for (i = 0; i < m; i++) { 267 for (j = 0; j < m; j++) Xi[i * m + j] = X[i * m + j] = (PetscScalar)(i == j); 268 speeds[i] = PETSC_MAX_REAL; /* Indicates invalid */ 269 } 270 PetscFunctionReturn(PETSC_SUCCESS); 271 } 272 273 static PetscErrorCode PhysicsCharacteristic_Shallow(void *vctx, PetscInt m, const PetscScalar *u, PetscScalar *X, PetscScalar *Xi, PetscReal *speeds) 274 { 275 ShallowCtx *phys = (ShallowCtx *)vctx; 276 PetscReal c; 277 PetscReal tol = 1e-6; 278 279 PetscFunctionBeginUser; 280 c = PetscSqrtScalar(u[0] * phys->gravity); 281 282 if (u[0] < tol) { /*Use conservative variables*/ 283 X[0 * 2 + 0] = 1; 284 X[0 * 2 + 1] = 0; 285 X[1 * 2 + 0] = 0; 286 X[1 * 2 + 1] = 1; 287 } else { 288 speeds[0] = u[1] / u[0] - c; 289 speeds[1] = u[1] / u[0] + c; 290 X[0 * 2 + 0] = 1; 291 X[0 * 2 + 1] = speeds[0]; 292 X[1 * 2 + 0] = 1; 293 X[1 * 2 + 1] = speeds[1]; 294 } 295 296 PetscCall(PetscArraycpy(Xi, X, 4)); 297 PetscCall(PetscKernel_A_gets_inverse_A_2(Xi, 0, PETSC_FALSE, NULL)); 298 PetscFunctionReturn(PETSC_SUCCESS); 299 } 300 301 static PetscErrorCode PhysicsSample_Shallow(void *vctx, PetscInt initial, FVBCType bctype, PetscReal xmin, PetscReal xmax, PetscReal t, PetscReal x, PetscReal *u) 302 { 303 PetscFunctionBeginUser; 304 PetscCheck(t <= 0, PETSC_COMM_SELF, PETSC_ERR_SUP, "Exact solutions not implemented for t > 0"); 305 switch (initial) { 306 case 0: 307 u[0] = (x < 0) ? 2 : 1; /* Standard Dam Break Problem */ 308 u[1] = (x < 0) ? 0 : 0; 309 break; 310 case 1: 311 u[0] = (x < 10) ? 1 : 0.1; /*The Next 5 problems are standard Riemann problem tests */ 312 u[1] = (x < 10) ? 2.5 : 0; 313 break; 314 case 2: 315 u[0] = (x < 25) ? 1 : 1; 316 u[1] = (x < 25) ? -5 : 5; 317 break; 318 case 3: 319 u[0] = (x < 20) ? 1 : 1e-12; 320 u[1] = (x < 20) ? 0 : 0; 321 break; 322 case 4: 323 u[0] = (x < 30) ? 1e-12 : 1; 324 u[1] = (x < 30) ? 0 : 0; 325 break; 326 case 5: 327 u[0] = (x < 25) ? 0.1 : 0.1; 328 u[1] = (x < 25) ? -0.3 : 0.3; 329 break; 330 case 6: 331 u[0] = 1 + 0.5 * PetscSinReal(2 * PETSC_PI * x); 332 u[1] = 1 * u[0]; 333 break; 334 case 7: 335 if (x < -0.1) { 336 u[0] = 1e-9; 337 u[1] = 0.0; 338 } else if (x < 0.1) { 339 u[0] = 1.0; 340 u[1] = 0.0; 341 } else { 342 u[0] = 1e-9; 343 u[1] = 0.0; 344 } 345 break; 346 case 8: 347 if (x < -0.1) { 348 u[0] = 2; 349 u[1] = 0.0; 350 } else if (x < 0.1) { 351 u[0] = 3.0; 352 u[1] = 0.0; 353 } else { 354 u[0] = 2; 355 u[1] = 0.0; 356 } 357 break; 358 default: 359 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "unknown initial condition"); 360 } 361 PetscFunctionReturn(PETSC_SUCCESS); 362 } 363 364 /* Implements inflow conditions for the given initial conditions. Which conditions are actually enforced depends on 365 on the results of PhysicsSetInflowType_Shallow. */ 366 static PetscErrorCode PhysicsInflow_Shallow(void *vctx, PetscReal t, PetscReal x, PetscReal *u) 367 { 368 FVCtx *ctx = (FVCtx *)vctx; 369 370 PetscFunctionBeginUser; 371 if (ctx->bctype == FVBC_INFLOW) { 372 switch (ctx->initial) { 373 case 0: 374 case 1: 375 case 2: 376 case 3: 377 case 4: 378 case 5: 379 u[0] = 0; 380 u[1] = 0.0; /* Left boundary conditions */ 381 u[2] = 0; 382 u[3] = 0.0; /* Right boundary conditions */ 383 break; 384 case 6: 385 u[0] = 0; 386 u[1] = 0.0; /* Left boundary conditions */ 387 u[2] = 0; 388 u[3] = 0.0; /* Right boundary conditions */ 389 break; 390 case 7: 391 u[0] = 0; 392 u[1] = 0.0; /* Left boundary conditions */ 393 u[2] = 0; 394 u[3] = 0.0; /* Right boundary conditions */ 395 break; 396 case 8: 397 u[0] = 0; 398 u[1] = 1.0; /* Left boundary conditions */ 399 u[2] = 0; 400 u[3] = -1.0; /* Right boundary conditions */ 401 break; 402 default: 403 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "unknown initial condition"); 404 } 405 } 406 PetscFunctionReturn(PETSC_SUCCESS); 407 } 408 409 /* Selects which boundary conditions are marked as inflow and which as outflow when FVBC_INFLOW is selected. */ 410 static PetscErrorCode PhysicsSetInflowType_Shallow(FVCtx *ctx) 411 { 412 PetscFunctionBeginUser; 413 switch (ctx->initial) { 414 case 0: 415 case 1: 416 case 2: 417 case 3: 418 case 4: 419 case 5: 420 case 6: 421 case 7: 422 case 8: /* Fix left and right momentum, height is outflow */ 423 ctx->physics2.bcinflowindex[0] = PETSC_FALSE; 424 ctx->physics2.bcinflowindex[1] = PETSC_TRUE; 425 ctx->physics2.bcinflowindex[2] = PETSC_FALSE; 426 ctx->physics2.bcinflowindex[3] = PETSC_TRUE; 427 break; 428 default: 429 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "unknown initial condition"); 430 } 431 PetscFunctionReturn(PETSC_SUCCESS); 432 } 433 434 static PetscErrorCode PhysicsCreate_Shallow(FVCtx *ctx) 435 { 436 ShallowCtx *user; 437 PetscFunctionList rlist = 0, rclist = 0; 438 char rname[256] = "rusanov", rcname[256] = "characteristic"; 439 440 PetscFunctionBeginUser; 441 PetscCall(PetscNew(&user)); 442 ctx->physics2.sample2 = PhysicsSample_Shallow; 443 ctx->physics2.inflow = PhysicsInflow_Shallow; 444 ctx->physics2.destroy = PhysicsDestroy_SimpleFree; 445 ctx->physics2.riemann2 = PhysicsRiemann_Shallow_Rusanov; 446 ctx->physics2.characteristic2 = PhysicsCharacteristic_Shallow; 447 ctx->physics2.user = user; 448 ctx->physics2.dof = 2; 449 450 PetscCall(PetscMalloc1(2 * (ctx->physics2.dof), &ctx->physics2.bcinflowindex)); 451 PetscCall(PhysicsSetInflowType_Shallow(ctx)); 452 453 PetscCall(PetscStrallocpy("height", &ctx->physics2.fieldname[0])); 454 PetscCall(PetscStrallocpy("momentum", &ctx->physics2.fieldname[1])); 455 456 user->gravity = 9.81; 457 458 PetscCall(RiemannListAdd_2WaySplit(&rlist, "exact", PhysicsRiemann_Shallow_Exact)); 459 PetscCall(RiemannListAdd_2WaySplit(&rlist, "rusanov", PhysicsRiemann_Shallow_Rusanov)); 460 PetscCall(ReconstructListAdd_2WaySplit(&rclist, "characteristic", PhysicsCharacteristic_Shallow)); 461 PetscCall(ReconstructListAdd_2WaySplit(&rclist, "conservative", PhysicsCharacteristic_Conservative)); 462 PetscOptionsBegin(ctx->comm, ctx->prefix, "Options for Shallow", ""); 463 PetscCall(PetscOptionsReal("-physics_shallow_gravity", "Gravity", "", user->gravity, &user->gravity, NULL)); 464 PetscCall(PetscOptionsFList("-physics_shallow_riemann", "Riemann solver", "", rlist, rname, rname, sizeof(rname), NULL)); 465 PetscCall(PetscOptionsFList("-physics_shallow_reconstruct", "Reconstruction", "", rclist, rcname, rcname, sizeof(rcname), NULL)); 466 PetscOptionsEnd(); 467 PetscCall(RiemannListFind_2WaySplit(rlist, rname, &ctx->physics2.riemann2)); 468 PetscCall(ReconstructListFind_2WaySplit(rclist, rcname, &ctx->physics2.characteristic2)); 469 PetscCall(PetscFunctionListDestroy(&rlist)); 470 PetscCall(PetscFunctionListDestroy(&rclist)); 471 PetscFunctionReturn(PETSC_SUCCESS); 472 } 473 474 PetscErrorCode FVSample_2WaySplit(FVCtx *ctx, DM da, PetscReal time, Vec U) 475 { 476 PetscScalar *u, *uj, xj, xi; 477 PetscInt i, j, k, dof, xs, xm, Mx; 478 const PetscInt N = 200; 479 PetscReal hs, hf; 480 481 PetscFunctionBeginUser; 482 PetscCheck(ctx->physics2.sample2, PETSC_COMM_SELF, PETSC_ERR_SUP, "Physics has not provided a sampling function"); 483 PetscCall(DMDAGetInfo(da, 0, &Mx, 0, 0, 0, 0, 0, &dof, 0, 0, 0, 0, 0)); 484 PetscCall(DMDAGetCorners(da, &xs, 0, 0, &xm, 0, 0)); 485 PetscCall(DMDAVecGetArray(da, U, &u)); 486 PetscCall(PetscMalloc1(dof, &uj)); 487 hs = (ctx->xmax - ctx->xmin) * 3.0 / 8.0 / ctx->sf; 488 hf = (ctx->xmax - ctx->xmin) / 4.0 / (ctx->fs - ctx->sf); 489 for (i = xs; i < xs + xm; i++) { 490 if (i < ctx->sf) { 491 xi = ctx->xmin + 0.5 * hs + i * hs; 492 /* Integrate over cell i using trapezoid rule with N points. */ 493 for (k = 0; k < dof; k++) u[i * dof + k] = 0; 494 for (j = 0; j < N + 1; j++) { 495 xj = xi + hs * (j - N / 2) / (PetscReal)N; 496 PetscCall((*ctx->physics2.sample2)(ctx->physics2.user, ctx->initial, ctx->bctype, ctx->xmin, ctx->xmax, time, xj, uj)); 497 for (k = 0; k < dof; k++) u[i * dof + k] += ((j == 0 || j == N) ? 0.5 : 1.0) * uj[k] / N; 498 } 499 } else if (i < ctx->fs) { 500 xi = ctx->xmin + ctx->sf * hs + 0.5 * hf + (i - ctx->sf) * hf; 501 /* Integrate over cell i using trapezoid rule with N points. */ 502 for (k = 0; k < dof; k++) u[i * dof + k] = 0; 503 for (j = 0; j < N + 1; j++) { 504 xj = xi + hf * (j - N / 2) / (PetscReal)N; 505 PetscCall((*ctx->physics2.sample2)(ctx->physics2.user, ctx->initial, ctx->bctype, ctx->xmin, ctx->xmax, time, xj, uj)); 506 for (k = 0; k < dof; k++) u[i * dof + k] += ((j == 0 || j == N) ? 0.5 : 1.0) * uj[k] / N; 507 } 508 } else { 509 xi = ctx->xmin + ctx->sf * hs + (ctx->fs - ctx->sf) * hf + 0.5 * hs + (i - ctx->fs) * hs; 510 /* Integrate over cell i using trapezoid rule with N points. */ 511 for (k = 0; k < dof; k++) u[i * dof + k] = 0; 512 for (j = 0; j < N + 1; j++) { 513 xj = xi + hs * (j - N / 2) / (PetscReal)N; 514 PetscCall((*ctx->physics2.sample2)(ctx->physics2.user, ctx->initial, ctx->bctype, ctx->xmin, ctx->xmax, time, xj, uj)); 515 for (k = 0; k < dof; k++) u[i * dof + k] += ((j == 0 || j == N) ? 0.5 : 1.0) * uj[k] / N; 516 } 517 } 518 } 519 PetscCall(DMDAVecRestoreArray(da, U, &u)); 520 PetscCall(PetscFree(uj)); 521 PetscFunctionReturn(PETSC_SUCCESS); 522 } 523 524 static PetscErrorCode SolutionErrorNorms_2WaySplit(FVCtx *ctx, DM da, PetscReal t, Vec X, PetscReal *nrm1) 525 { 526 Vec Y; 527 PetscInt i, Mx; 528 const PetscScalar *ptr_X, *ptr_Y; 529 PetscReal hs, hf; 530 531 PetscFunctionBeginUser; 532 PetscCall(VecGetSize(X, &Mx)); 533 PetscCall(VecDuplicate(X, &Y)); 534 PetscCall(FVSample_2WaySplit(ctx, da, t, Y)); 535 hs = (ctx->xmax - ctx->xmin) * 3.0 / 8.0 / ctx->sf; 536 hf = (ctx->xmax - ctx->xmin) / 4.0 / (ctx->fs - ctx->sf); 537 PetscCall(VecGetArrayRead(X, &ptr_X)); 538 PetscCall(VecGetArrayRead(Y, &ptr_Y)); 539 for (i = 0; i < Mx; i++) { 540 if (i < ctx->sf || i > ctx->fs - 1) *nrm1 += hs * PetscAbs(ptr_X[i] - ptr_Y[i]); 541 else *nrm1 += hf * PetscAbs(ptr_X[i] - ptr_Y[i]); 542 } 543 PetscCall(VecRestoreArrayRead(X, &ptr_X)); 544 PetscCall(VecRestoreArrayRead(Y, &ptr_Y)); 545 PetscCall(VecDestroy(&Y)); 546 PetscFunctionReturn(PETSC_SUCCESS); 547 } 548 549 PetscErrorCode FVRHSFunction_2WaySplit(TS ts, PetscReal time, Vec X, Vec F, void *vctx) 550 { 551 FVCtx *ctx = (FVCtx *)vctx; 552 PetscInt i, j, k, Mx, dof, xs, xm, sf = ctx->sf, fs = ctx->fs; 553 PetscReal hxf, hxs, cfl_idt = 0; 554 PetscScalar *x, *f, *slope; 555 Vec Xloc; 556 DM da; 557 558 PetscFunctionBeginUser; 559 PetscCall(TSGetDM(ts, &da)); 560 PetscCall(DMGetLocalVector(da, &Xloc)); /* Xloc contains ghost points */ 561 PetscCall(DMDAGetInfo(da, 0, &Mx, 0, 0, 0, 0, 0, &dof, 0, 0, 0, 0, 0)); /* Mx is the number of center points */ 562 hxs = (ctx->xmax - ctx->xmin) * 3.0 / 8.0 / ctx->sf; 563 hxf = (ctx->xmax - ctx->xmin) / 4.0 / (ctx->fs - ctx->sf); 564 PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, Xloc)); /* X is solution vector which does not contain ghost points */ 565 PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, Xloc)); 566 567 PetscCall(VecZeroEntries(F)); /* F is the right-hand side function corresponds to center points */ 568 569 PetscCall(DMDAVecGetArray(da, Xloc, &x)); 570 PetscCall(DMDAVecGetArray(da, F, &f)); 571 PetscCall(DMDAGetArray(da, PETSC_TRUE, &slope)); /* contains ghost points */ 572 PetscCall(DMDAGetCorners(da, &xs, 0, 0, &xm, 0, 0)); 573 574 if (ctx->bctype == FVBC_OUTFLOW) { 575 for (i = xs - 2; i < 0; i++) { 576 for (j = 0; j < dof; j++) x[i * dof + j] = x[j]; 577 } 578 for (i = Mx; i < xs + xm + 2; i++) { 579 for (j = 0; j < dof; j++) x[i * dof + j] = x[(xs + xm - 1) * dof + j]; 580 } 581 } 582 583 if (ctx->bctype == FVBC_INFLOW) { 584 /* See LeVeque, R. (2002). Finite Volume Methods for Hyperbolic Problems. doi:10.1017/CBO9780511791253 585 pages 137-138 for the scheme. */ 586 if (xs == 0) { /* Left Boundary */ 587 PetscCall(ctx->physics2.inflow(ctx, time, ctx->xmin, ctx->ub)); 588 for (j = 0; j < dof; j++) { 589 if (ctx->physics2.bcinflowindex[j]) { 590 for (i = -2; i < 0; i++) x[i * dof + j] = 2.0 * ctx->ub[j] - x[-(i + 1) * dof + j]; 591 } else { 592 for (i = -2; i < 0; i++) x[i * dof + j] = x[j]; /* Outflow */ 593 } 594 } 595 } 596 if (xs + xm == Mx) { /* Right Boundary */ 597 PetscCall(ctx->physics2.inflow(ctx, time, ctx->xmax, ctx->ub)); 598 for (j = 0; j < dof; j++) { 599 if (ctx->physics2.bcinflowindex[dof + j]) { 600 for (i = Mx; i < Mx + 2; i++) x[i * dof + j] = 2.0 * ctx->ub[dof + j] - x[(2 * Mx - (i + 1)) * dof + j]; 601 } else { 602 for (i = Mx; i < Mx + 2; i++) x[i * dof + j] = x[(Mx - 1) * dof + j]; /* Outflow */ 603 } 604 } 605 } 606 } 607 608 for (i = xs - 1; i < xs + xm + 1; i++) { 609 struct _LimitInfo info; 610 PetscScalar *cjmpL, *cjmpR; 611 /* Determine the right eigenvectors R, where A = R \Lambda R^{-1} */ 612 PetscCall((*ctx->physics2.characteristic2)(ctx->physics2.user, dof, &x[i * dof], ctx->R, ctx->Rinv, ctx->speeds)); 613 /* Evaluate jumps across interfaces (i-1, i) and (i, i+1), put in characteristic basis */ 614 PetscCall(PetscArrayzero(ctx->cjmpLR, 2 * dof)); 615 cjmpL = &ctx->cjmpLR[0]; 616 cjmpR = &ctx->cjmpLR[dof]; 617 for (j = 0; j < dof; j++) { 618 PetscScalar jmpL, jmpR; 619 jmpL = x[(i + 0) * dof + j] - x[(i - 1) * dof + j]; 620 jmpR = x[(i + 1) * dof + j] - x[(i + 0) * dof + j]; 621 for (k = 0; k < dof; k++) { 622 cjmpL[k] += ctx->Rinv[k + j * dof] * jmpL; 623 cjmpR[k] += ctx->Rinv[k + j * dof] * jmpR; 624 } 625 } 626 /* Apply limiter to the left and right characteristic jumps */ 627 info.m = dof; 628 info.hxs = hxs; 629 info.hxf = hxf; 630 (*ctx->limit2)(&info, cjmpL, cjmpR, ctx->sf, ctx->fs, i, ctx->cslope); 631 for (j = 0; j < dof; j++) { 632 PetscScalar tmp = 0; 633 for (k = 0; k < dof; k++) tmp += ctx->R[j + k * dof] * ctx->cslope[k]; 634 slope[i * dof + j] = tmp; 635 } 636 } 637 638 for (i = xs; i < xs + xm + 1; i++) { 639 PetscReal maxspeed; 640 PetscScalar *uL, *uR; 641 uL = &ctx->uLR[0]; 642 uR = &ctx->uLR[dof]; 643 if (i < sf) { /* slow region */ 644 for (j = 0; j < dof; j++) { 645 uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2; 646 uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2; 647 } 648 PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed)); 649 if (i > xs) { 650 for (j = 0; j < dof; j++) f[(i - 1) * dof + j] -= ctx->flux[j] / hxs; 651 } 652 if (i < xs + xm) { 653 for (j = 0; j < dof; j++) f[i * dof + j] += ctx->flux[j] / hxs; 654 } 655 } else if (i == sf) { /* interface between the slow region and the fast region */ 656 for (j = 0; j < dof; j++) { 657 uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2; 658 uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxf / 2; 659 } 660 PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed)); 661 if (i > xs) { 662 for (j = 0; j < dof; j++) f[(i - 1) * dof + j] -= ctx->flux[j] / hxs; 663 } 664 if (i < xs + xm) { 665 for (j = 0; j < dof; j++) f[i * dof + j] += ctx->flux[j] / hxf; 666 } 667 } else if (i > sf && i < fs) { /* fast region */ 668 for (j = 0; j < dof; j++) { 669 uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxf / 2; 670 uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxf / 2; 671 } 672 PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed)); 673 if (i > xs) { 674 for (j = 0; j < dof; j++) f[(i - 1) * dof + j] -= ctx->flux[j] / hxf; 675 } 676 if (i < xs + xm) { 677 for (j = 0; j < dof; j++) f[i * dof + j] += ctx->flux[j] / hxf; 678 } 679 } else if (i == fs) { /* interface between the fast region and the slow region */ 680 for (j = 0; j < dof; j++) { 681 uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxf / 2; 682 uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2; 683 } 684 PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed)); 685 if (i > xs) { 686 for (j = 0; j < dof; j++) f[(i - 1) * dof + j] -= ctx->flux[j] / hxf; 687 } 688 if (i < xs + xm) { 689 for (j = 0; j < dof; j++) f[i * dof + j] += ctx->flux[j] / hxs; 690 } 691 } else { /* slow region */ 692 for (j = 0; j < dof; j++) { 693 uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2; 694 uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2; 695 } 696 PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed)); 697 cfl_idt = PetscMax(cfl_idt, PetscAbsScalar(maxspeed / hxs)); /* Max allowable value of 1/Delta t */ 698 if (i > xs) { 699 for (j = 0; j < dof; j++) f[(i - 1) * dof + j] -= ctx->flux[j] / hxs; 700 } 701 if (i < xs + xm) { 702 for (j = 0; j < dof; j++) f[i * dof + j] += ctx->flux[j] / hxs; 703 } 704 } 705 } 706 PetscCall(DMDAVecRestoreArray(da, Xloc, &x)); 707 PetscCall(DMDAVecRestoreArray(da, F, &f)); 708 PetscCall(DMDARestoreArray(da, PETSC_TRUE, &slope)); 709 PetscCall(DMRestoreLocalVector(da, &Xloc)); 710 PetscCall(MPIU_Allreduce(&cfl_idt, &ctx->cfl_idt, 1, MPIU_SCALAR, MPIU_MAX, PetscObjectComm((PetscObject)da))); 711 if (0) { 712 /* We need to a way to inform the TS of a CFL constraint, this is a debugging fragment */ 713 PetscReal dt, tnow; 714 PetscCall(TSGetTimeStep(ts, &dt)); 715 PetscCall(TSGetTime(ts, &tnow)); 716 if (dt > 0.5 / ctx->cfl_idt) { 717 if (1) { 718 PetscCall(PetscPrintf(ctx->comm, "Stability constraint exceeded at t=%g, dt %g > %g\n", (double)tnow, (double)dt, (double)(0.5 / ctx->cfl_idt))); 719 } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Stability constraint exceeded, %g > %g", (double)dt, (double)(ctx->cfl / ctx->cfl_idt)); 720 } 721 } 722 PetscFunctionReturn(PETSC_SUCCESS); 723 } 724 725 /* --------------------------------- Finite Volume Solver for slow components ----------------------------------- */ 726 PetscErrorCode FVRHSFunctionslow_2WaySplit(TS ts, PetscReal time, Vec X, Vec F, void *vctx) 727 { 728 FVCtx *ctx = (FVCtx *)vctx; 729 PetscInt i, j, k, Mx, dof, xs, xm, islow = 0, sf = ctx->sf, fs = ctx->fs, lsbwidth = ctx->lsbwidth, rsbwidth = ctx->rsbwidth; 730 PetscReal hxs, hxf, cfl_idt = 0; 731 PetscScalar *x, *f, *slope; 732 Vec Xloc; 733 DM da; 734 735 PetscFunctionBeginUser; 736 PetscCall(TSGetDM(ts, &da)); 737 PetscCall(DMGetLocalVector(da, &Xloc)); 738 PetscCall(DMDAGetInfo(da, 0, &Mx, 0, 0, 0, 0, 0, &dof, 0, 0, 0, 0, 0)); 739 hxs = (ctx->xmax - ctx->xmin) * 3.0 / 8.0 / ctx->sf; 740 hxf = (ctx->xmax - ctx->xmin) / 4.0 / (ctx->fs - ctx->sf); 741 PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, Xloc)); 742 PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, Xloc)); 743 PetscCall(VecZeroEntries(F)); 744 PetscCall(DMDAVecGetArray(da, Xloc, &x)); 745 PetscCall(VecGetArray(F, &f)); 746 PetscCall(DMDAGetArray(da, PETSC_TRUE, &slope)); 747 PetscCall(DMDAGetCorners(da, &xs, 0, 0, &xm, 0, 0)); 748 749 if (ctx->bctype == FVBC_OUTFLOW) { 750 for (i = xs - 2; i < 0; i++) { 751 for (j = 0; j < dof; j++) x[i * dof + j] = x[j]; 752 } 753 for (i = Mx; i < xs + xm + 2; i++) { 754 for (j = 0; j < dof; j++) x[i * dof + j] = x[(xs + xm - 1) * dof + j]; 755 } 756 } 757 if (ctx->bctype == FVBC_INFLOW) { 758 /* See LeVeque, R. (2002). Finite Volume Methods for Hyperbolic Problems. doi:10.1017/CBO9780511791253 759 pages 137-138 for the scheme. */ 760 if (xs == 0) { /* Left Boundary */ 761 PetscCall(ctx->physics2.inflow(ctx, time, ctx->xmin, ctx->ub)); 762 for (j = 0; j < dof; j++) { 763 if (ctx->physics2.bcinflowindex[j] == PETSC_TRUE) { 764 for (i = -2; i < 0; i++) x[i * dof + j] = 2.0 * ctx->ub[j] - x[-(i + 1) * dof + j]; 765 } else { 766 for (i = -2; i < 0; i++) x[i * dof + j] = x[j]; /* Outflow */ 767 } 768 } 769 } 770 if (xs + xm == Mx) { /* Right Boundary */ 771 PetscCall(ctx->physics2.inflow(ctx, time, ctx->xmax, ctx->ub)); 772 for (j = 0; j < dof; j++) { 773 if (ctx->physics2.bcinflowindex[dof + j] == PETSC_TRUE) { 774 for (i = Mx; i < Mx + 2; i++) x[i * dof + j] = 2.0 * ctx->ub[dof + j] - x[(2 * Mx - (i + 1)) * dof + j]; 775 } else { 776 for (i = Mx; i < Mx + 2; i++) x[i * dof + j] = x[(Mx - 1) * dof + j]; /* Outflow */ 777 } 778 } 779 } 780 } 781 for (i = xs - 1; i < xs + xm + 1; i++) { 782 struct _LimitInfo info; 783 PetscScalar *cjmpL, *cjmpR; 784 if (i < sf - lsbwidth + 1 || i > fs + rsbwidth - 2) { /* slow components and the first and last fast components */ 785 /* Determine the right eigenvectors R, where A = R \Lambda R^{-1} */ 786 PetscCall((*ctx->physics2.characteristic2)(ctx->physics2.user, dof, &x[i * dof], ctx->R, ctx->Rinv, ctx->speeds)); 787 /* Evaluate jumps across interfaces (i-1, i) and (i, i+1), put in characteristic basis */ 788 PetscCall(PetscArrayzero(ctx->cjmpLR, 2 * dof)); 789 cjmpL = &ctx->cjmpLR[0]; 790 cjmpR = &ctx->cjmpLR[dof]; 791 for (j = 0; j < dof; j++) { 792 PetscScalar jmpL, jmpR; 793 jmpL = x[(i + 0) * dof + j] - x[(i - 1) * dof + j]; 794 jmpR = x[(i + 1) * dof + j] - x[(i + 0) * dof + j]; 795 for (k = 0; k < dof; k++) { 796 cjmpL[k] += ctx->Rinv[k + j * dof] * jmpL; 797 cjmpR[k] += ctx->Rinv[k + j * dof] * jmpR; 798 } 799 } 800 /* Apply limiter to the left and right characteristic jumps */ 801 info.m = dof; 802 info.hxs = hxs; 803 info.hxf = hxf; 804 (*ctx->limit2)(&info, cjmpL, cjmpR, ctx->sf, ctx->fs, i, ctx->cslope); 805 for (j = 0; j < dof; j++) { 806 PetscScalar tmp = 0; 807 for (k = 0; k < dof; k++) tmp += ctx->R[j + k * dof] * ctx->cslope[k]; 808 slope[i * dof + j] = tmp; 809 } 810 } 811 } 812 813 for (i = xs; i < xs + xm + 1; i++) { 814 PetscReal maxspeed; 815 PetscScalar *uL, *uR; 816 uL = &ctx->uLR[0]; 817 uR = &ctx->uLR[dof]; 818 if (i < sf - lsbwidth) { /* slow region */ 819 for (j = 0; j < dof; j++) { 820 uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2; 821 uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2; 822 } 823 PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed)); 824 cfl_idt = PetscMax(cfl_idt, PetscAbsScalar(maxspeed / hxs)); /* Max allowable value of 1/Delta t */ 825 if (i > xs) { 826 for (j = 0; j < dof; j++) f[(islow - 1) * dof + j] -= ctx->flux[j] / hxs; 827 } 828 if (i < xs + xm) { 829 for (j = 0; j < dof; j++) f[islow * dof + j] += ctx->flux[j] / hxs; 830 islow++; 831 } 832 } 833 if (i == sf - lsbwidth) { /* interface between the slow region and the fast region */ 834 for (j = 0; j < dof; j++) { 835 uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2; 836 uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2; 837 } 838 PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed)); 839 if (i > xs) { 840 for (j = 0; j < dof; j++) f[(islow - 1) * dof + j] -= ctx->flux[j] / hxs; 841 } 842 } 843 if (i == fs + rsbwidth) { /* slow region */ 844 for (j = 0; j < dof; j++) { 845 uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2; 846 uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2; 847 } 848 PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed)); 849 if (i < xs + xm) { 850 for (j = 0; j < dof; j++) f[islow * dof + j] += ctx->flux[j] / hxs; 851 islow++; 852 } 853 } 854 if (i > fs + rsbwidth) { /* slow region */ 855 for (j = 0; j < dof; j++) { 856 uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2; 857 uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2; 858 } 859 PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed)); 860 if (i > xs) { 861 for (j = 0; j < dof; j++) f[(islow - 1) * dof + j] -= ctx->flux[j] / hxs; 862 } 863 if (i < xs + xm) { 864 for (j = 0; j < dof; j++) f[islow * dof + j] += ctx->flux[j] / hxs; 865 islow++; 866 } 867 } 868 } 869 PetscCall(DMDAVecRestoreArray(da, Xloc, &x)); 870 PetscCall(VecRestoreArray(F, &f)); 871 PetscCall(DMDARestoreArray(da, PETSC_TRUE, &slope)); 872 PetscCall(DMRestoreLocalVector(da, &Xloc)); 873 PetscCall(MPIU_Allreduce(&cfl_idt, &ctx->cfl_idt, 1, MPIU_SCALAR, MPIU_MAX, PetscObjectComm((PetscObject)da))); 874 PetscFunctionReturn(PETSC_SUCCESS); 875 } 876 877 PetscErrorCode FVRHSFunctionslowbuffer_2WaySplit(TS ts, PetscReal time, Vec X, Vec F, void *vctx) 878 { 879 FVCtx *ctx = (FVCtx *)vctx; 880 PetscInt i, j, k, Mx, dof, xs, xm, islow = 0, sf = ctx->sf, fs = ctx->fs, lsbwidth = ctx->lsbwidth, rsbwidth = ctx->rsbwidth; 881 PetscReal hxs, hxf; 882 PetscScalar *x, *f, *slope; 883 Vec Xloc; 884 DM da; 885 886 PetscFunctionBeginUser; 887 PetscCall(TSGetDM(ts, &da)); 888 PetscCall(DMGetLocalVector(da, &Xloc)); 889 PetscCall(DMDAGetInfo(da, 0, &Mx, 0, 0, 0, 0, 0, &dof, 0, 0, 0, 0, 0)); 890 hxs = (ctx->xmax - ctx->xmin) * 3.0 / 8.0 / ctx->sf; 891 hxf = (ctx->xmax - ctx->xmin) / 4.0 / (ctx->fs - ctx->sf); 892 PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, Xloc)); 893 PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, Xloc)); 894 PetscCall(VecZeroEntries(F)); 895 PetscCall(DMDAVecGetArray(da, Xloc, &x)); 896 PetscCall(VecGetArray(F, &f)); 897 PetscCall(DMDAGetArray(da, PETSC_TRUE, &slope)); 898 PetscCall(DMDAGetCorners(da, &xs, 0, 0, &xm, 0, 0)); 899 900 if (ctx->bctype == FVBC_OUTFLOW) { 901 for (i = xs - 2; i < 0; i++) { 902 for (j = 0; j < dof; j++) x[i * dof + j] = x[j]; 903 } 904 for (i = Mx; i < xs + xm + 2; i++) { 905 for (j = 0; j < dof; j++) x[i * dof + j] = x[(xs + xm - 1) * dof + j]; 906 } 907 } 908 if (ctx->bctype == FVBC_INFLOW) { 909 /* See LeVeque, R. (2002). Finite Volume Methods for Hyperbolic Problems. doi:10.1017/CBO9780511791253 910 pages 137-138 for the scheme. */ 911 if (xs == 0) { /* Left Boundary */ 912 PetscCall(ctx->physics2.inflow(ctx, time, ctx->xmin, ctx->ub)); 913 for (j = 0; j < dof; j++) { 914 if (ctx->physics2.bcinflowindex[j] == PETSC_TRUE) { 915 for (i = -2; i < 0; i++) x[i * dof + j] = 2.0 * ctx->ub[j] - x[-(i + 1) * dof + j]; 916 } else { 917 for (i = -2; i < 0; i++) x[i * dof + j] = x[j]; /* Outflow */ 918 } 919 } 920 } 921 if (xs + xm == Mx) { /* Right Boundary */ 922 PetscCall(ctx->physics2.inflow(ctx, time, ctx->xmax, ctx->ub)); 923 for (j = 0; j < dof; j++) { 924 if (ctx->physics2.bcinflowindex[dof + j] == PETSC_TRUE) { 925 for (i = Mx; i < Mx + 2; i++) x[i * dof + j] = 2.0 * ctx->ub[dof + j] - x[(2 * Mx - (i + 1)) * dof + j]; 926 } else { 927 for (i = Mx; i < Mx + 2; i++) x[i * dof + j] = x[(Mx - 1) * dof + j]; /* Outflow */ 928 } 929 } 930 } 931 } 932 for (i = xs - 1; i < xs + xm + 1; i++) { 933 struct _LimitInfo info; 934 PetscScalar *cjmpL, *cjmpR; 935 if ((i > sf - lsbwidth - 2 && i < sf + 1) || (i > fs - 2 && i < fs + rsbwidth + 1)) { 936 /* Determine the right eigenvectors R, where A = R \Lambda R^{-1} */ 937 PetscCall((*ctx->physics2.characteristic2)(ctx->physics2.user, dof, &x[i * dof], ctx->R, ctx->Rinv, ctx->speeds)); 938 /* Evaluate jumps across interfaces (i-1, i) and (i, i+1), put in characteristic basis */ 939 PetscCall(PetscArrayzero(ctx->cjmpLR, 2 * dof)); 940 cjmpL = &ctx->cjmpLR[0]; 941 cjmpR = &ctx->cjmpLR[dof]; 942 for (j = 0; j < dof; j++) { 943 PetscScalar jmpL, jmpR; 944 jmpL = x[(i + 0) * dof + j] - x[(i - 1) * dof + j]; 945 jmpR = x[(i + 1) * dof + j] - x[(i + 0) * dof + j]; 946 for (k = 0; k < dof; k++) { 947 cjmpL[k] += ctx->Rinv[k + j * dof] * jmpL; 948 cjmpR[k] += ctx->Rinv[k + j * dof] * jmpR; 949 } 950 } 951 /* Apply limiter to the left and right characteristic jumps */ 952 info.m = dof; 953 info.hxs = hxs; 954 info.hxf = hxf; 955 (*ctx->limit2)(&info, cjmpL, cjmpR, ctx->sf, ctx->fs, i, ctx->cslope); 956 for (j = 0; j < dof; j++) { 957 PetscScalar tmp = 0; 958 for (k = 0; k < dof; k++) tmp += ctx->R[j + k * dof] * ctx->cslope[k]; 959 slope[i * dof + j] = tmp; 960 } 961 } 962 } 963 964 for (i = xs; i < xs + xm + 1; i++) { 965 PetscReal maxspeed; 966 PetscScalar *uL, *uR; 967 uL = &ctx->uLR[0]; 968 uR = &ctx->uLR[dof]; 969 if (i == sf - lsbwidth) { 970 for (j = 0; j < dof; j++) { 971 uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2; 972 uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2; 973 } 974 PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed)); 975 if (i < xs + xm) { 976 for (j = 0; j < dof; j++) f[islow * dof + j] += ctx->flux[j] / hxs; 977 islow++; 978 } 979 } 980 if (i > sf - lsbwidth && i < sf) { 981 for (j = 0; j < dof; j++) { 982 uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2; 983 uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2; 984 } 985 PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed)); 986 if (i > xs) { 987 for (j = 0; j < dof; j++) f[(islow - 1) * dof + j] -= ctx->flux[j] / hxs; 988 } 989 if (i < xs + xm) { 990 for (j = 0; j < dof; j++) f[islow * dof + j] += ctx->flux[j] / hxs; 991 islow++; 992 } 993 } 994 if (i == sf) { /* interface between the slow region and the fast region */ 995 for (j = 0; j < dof; j++) { 996 uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2; 997 uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxf / 2; 998 } 999 PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed)); 1000 if (i > xs) { 1001 for (j = 0; j < dof; j++) f[(islow - 1) * dof + j] -= ctx->flux[j] / hxs; 1002 } 1003 } 1004 if (i == fs) { /* interface between the fast region and the slow region */ 1005 for (j = 0; j < dof; j++) { 1006 uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxf / 2; 1007 uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2; 1008 } 1009 PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed)); 1010 if (i < xs + xm) { 1011 for (j = 0; j < dof; j++) f[islow * dof + j] += ctx->flux[j] / hxs; 1012 islow++; 1013 } 1014 } 1015 if (i > fs && i < fs + rsbwidth) { 1016 for (j = 0; j < dof; j++) { 1017 uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2; 1018 uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2; 1019 } 1020 PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed)); 1021 if (i > xs) { 1022 for (j = 0; j < dof; j++) f[(islow - 1) * dof + j] -= ctx->flux[j] / hxs; 1023 } 1024 if (i < xs + xm) { 1025 for (j = 0; j < dof; j++) f[islow * dof + j] += ctx->flux[j] / hxs; 1026 islow++; 1027 } 1028 } 1029 if (i == fs + rsbwidth) { 1030 for (j = 0; j < dof; j++) { 1031 uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2; 1032 uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2; 1033 } 1034 PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed)); 1035 if (i > xs) { 1036 for (j = 0; j < dof; j++) f[(islow - 1) * dof + j] -= ctx->flux[j] / hxs; 1037 } 1038 } 1039 } 1040 PetscCall(DMDAVecRestoreArray(da, Xloc, &x)); 1041 PetscCall(VecRestoreArray(F, &f)); 1042 PetscCall(DMDARestoreArray(da, PETSC_TRUE, &slope)); 1043 PetscCall(DMRestoreLocalVector(da, &Xloc)); 1044 PetscFunctionReturn(PETSC_SUCCESS); 1045 } 1046 1047 /* --------------------------------- Finite Volume Solver for fast parts ----------------------------------- */ 1048 PetscErrorCode FVRHSFunctionfast_2WaySplit(TS ts, PetscReal time, Vec X, Vec F, void *vctx) 1049 { 1050 FVCtx *ctx = (FVCtx *)vctx; 1051 PetscInt i, j, k, Mx, dof, xs, xm, ifast = 0, sf = ctx->sf, fs = ctx->fs; 1052 PetscReal hxs, hxf; 1053 PetscScalar *x, *f, *slope; 1054 Vec Xloc; 1055 DM da; 1056 1057 PetscFunctionBeginUser; 1058 PetscCall(TSGetDM(ts, &da)); 1059 PetscCall(DMGetLocalVector(da, &Xloc)); 1060 PetscCall(DMDAGetInfo(da, 0, &Mx, 0, 0, 0, 0, 0, &dof, 0, 0, 0, 0, 0)); 1061 hxs = (ctx->xmax - ctx->xmin) * 3.0 / 8.0 / ctx->sf; 1062 hxf = (ctx->xmax - ctx->xmin) / 4.0 / (ctx->fs - ctx->sf); 1063 PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, Xloc)); 1064 PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, Xloc)); 1065 PetscCall(VecZeroEntries(F)); 1066 PetscCall(DMDAVecGetArray(da, Xloc, &x)); 1067 PetscCall(VecGetArray(F, &f)); 1068 PetscCall(DMDAGetArray(da, PETSC_TRUE, &slope)); 1069 PetscCall(DMDAGetCorners(da, &xs, 0, 0, &xm, 0, 0)); 1070 1071 if (ctx->bctype == FVBC_OUTFLOW) { 1072 for (i = xs - 2; i < 0; i++) { 1073 for (j = 0; j < dof; j++) x[i * dof + j] = x[j]; 1074 } 1075 for (i = Mx; i < xs + xm + 2; i++) { 1076 for (j = 0; j < dof; j++) x[i * dof + j] = x[(xs + xm - 1) * dof + j]; 1077 } 1078 } 1079 if (ctx->bctype == FVBC_INFLOW) { 1080 /* See LeVeque, R. (2002). Finite Volume Methods for Hyperbolic Problems. doi:10.1017/CBO9780511791253 1081 pages 137-138 for the scheme. */ 1082 if (xs == 0) { /* Left Boundary */ 1083 PetscCall(ctx->physics2.inflow(ctx, time, ctx->xmin, ctx->ub)); 1084 for (j = 0; j < dof; j++) { 1085 if (ctx->physics2.bcinflowindex[j] == PETSC_TRUE) { 1086 for (i = -2; i < 0; i++) x[i * dof + j] = 2.0 * ctx->ub[j] - x[-(i + 1) * dof + j]; 1087 } else { 1088 for (i = -2; i < 0; i++) x[i * dof + j] = x[j]; /* Outflow */ 1089 } 1090 } 1091 } 1092 if (xs + xm == Mx) { /* Right Boundary */ 1093 PetscCall(ctx->physics2.inflow(ctx, time, ctx->xmax, ctx->ub)); 1094 for (j = 0; j < dof; j++) { 1095 if (ctx->physics2.bcinflowindex[dof + j] == PETSC_TRUE) { 1096 for (i = Mx; i < Mx + 2; i++) x[i * dof + j] = 2.0 * ctx->ub[dof + j] - x[(2 * Mx - (i + 1)) * dof + j]; 1097 } else { 1098 for (i = Mx; i < Mx + 2; i++) x[i * dof + j] = x[(Mx - 1) * dof + j]; /* Outflow */ 1099 } 1100 } 1101 } 1102 } 1103 for (i = xs - 1; i < xs + xm + 1; i++) { 1104 struct _LimitInfo info; 1105 PetscScalar *cjmpL, *cjmpR; 1106 if (i > sf - 2 && i < fs + 1) { 1107 PetscCall((*ctx->physics2.characteristic2)(ctx->physics2.user, dof, &x[i * dof], ctx->R, ctx->Rinv, ctx->speeds)); 1108 PetscCall(PetscArrayzero(ctx->cjmpLR, 2 * dof)); 1109 cjmpL = &ctx->cjmpLR[0]; 1110 cjmpR = &ctx->cjmpLR[dof]; 1111 for (j = 0; j < dof; j++) { 1112 PetscScalar jmpL, jmpR; 1113 jmpL = x[(i + 0) * dof + j] - x[(i - 1) * dof + j]; 1114 jmpR = x[(i + 1) * dof + j] - x[(i + 0) * dof + j]; 1115 for (k = 0; k < dof; k++) { 1116 cjmpL[k] += ctx->Rinv[k + j * dof] * jmpL; 1117 cjmpR[k] += ctx->Rinv[k + j * dof] * jmpR; 1118 } 1119 } 1120 /* Apply limiter to the left and right characteristic jumps */ 1121 info.m = dof; 1122 info.hxs = hxs; 1123 info.hxf = hxf; 1124 (*ctx->limit2)(&info, cjmpL, cjmpR, ctx->sf, ctx->fs, i, ctx->cslope); 1125 for (j = 0; j < dof; j++) { 1126 PetscScalar tmp = 0; 1127 for (k = 0; k < dof; k++) tmp += ctx->R[j + k * dof] * ctx->cslope[k]; 1128 slope[i * dof + j] = tmp; 1129 } 1130 } 1131 } 1132 1133 for (i = xs; i < xs + xm + 1; i++) { 1134 PetscReal maxspeed; 1135 PetscScalar *uL, *uR; 1136 uL = &ctx->uLR[0]; 1137 uR = &ctx->uLR[dof]; 1138 if (i == sf) { /* interface between the slow region and the fast region */ 1139 for (j = 0; j < dof; j++) { 1140 uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2; 1141 uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxf / 2; 1142 } 1143 PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed)); 1144 if (i < xs + xm) { 1145 for (j = 0; j < dof; j++) f[ifast * dof + j] += ctx->flux[j] / hxf; 1146 ifast++; 1147 } 1148 } 1149 if (i > sf && i < fs) { /* fast region */ 1150 for (j = 0; j < dof; j++) { 1151 uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxf / 2; 1152 uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxf / 2; 1153 } 1154 PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed)); 1155 if (i > xs) { 1156 for (j = 0; j < dof; j++) f[(ifast - 1) * dof + j] -= ctx->flux[j] / hxf; 1157 } 1158 if (i < xs + xm) { 1159 for (j = 0; j < dof; j++) f[ifast * dof + j] += ctx->flux[j] / hxf; 1160 ifast++; 1161 } 1162 } 1163 if (i == fs) { /* interface between the fast region and the slow region */ 1164 for (j = 0; j < dof; j++) { 1165 uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxf / 2; 1166 uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2; 1167 } 1168 PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed)); 1169 if (i > xs) { 1170 for (j = 0; j < dof; j++) f[(ifast - 1) * dof + j] -= ctx->flux[j] / hxf; 1171 } 1172 } 1173 } 1174 PetscCall(DMDAVecRestoreArray(da, Xloc, &x)); 1175 PetscCall(VecRestoreArray(F, &f)); 1176 PetscCall(DMDARestoreArray(da, PETSC_TRUE, &slope)); 1177 PetscCall(DMRestoreLocalVector(da, &Xloc)); 1178 PetscFunctionReturn(PETSC_SUCCESS); 1179 } 1180 1181 int main(int argc, char *argv[]) 1182 { 1183 char lname[256] = "mc", physname[256] = "advect", final_fname[256] = "solution.m"; 1184 PetscFunctionList limiters = 0, physics = 0; 1185 MPI_Comm comm; 1186 TS ts; 1187 DM da; 1188 Vec X, X0, R; 1189 FVCtx ctx; 1190 PetscInt i, k, dof, xs, xm, Mx, draw = 0, count_slow, count_fast, islow = 0, ifast = 0, islowbuffer = 0, *index_slow, *index_fast, *index_slowbuffer; 1191 PetscBool view_final = PETSC_FALSE; 1192 PetscReal ptime, maxtime; 1193 1194 PetscFunctionBeginUser; 1195 PetscCall(PetscInitialize(&argc, &argv, 0, help)); 1196 comm = PETSC_COMM_WORLD; 1197 PetscCall(PetscMemzero(&ctx, sizeof(ctx))); 1198 1199 /* Register limiters to be available on the command line */ 1200 PetscCall(PetscFunctionListAdd(&limiters, "upwind", Limit2_Upwind)); 1201 PetscCall(PetscFunctionListAdd(&limiters, "lax-wendroff", Limit2_LaxWendroff)); 1202 PetscCall(PetscFunctionListAdd(&limiters, "beam-warming", Limit2_BeamWarming)); 1203 PetscCall(PetscFunctionListAdd(&limiters, "fromm", Limit2_Fromm)); 1204 PetscCall(PetscFunctionListAdd(&limiters, "minmod", Limit2_Minmod)); 1205 PetscCall(PetscFunctionListAdd(&limiters, "superbee", Limit2_Superbee)); 1206 PetscCall(PetscFunctionListAdd(&limiters, "mc", Limit2_MC)); 1207 PetscCall(PetscFunctionListAdd(&limiters, "koren3", Limit2_Koren3)); 1208 1209 /* Register physical models to be available on the command line */ 1210 PetscCall(PetscFunctionListAdd(&physics, "shallow", PhysicsCreate_Shallow)); 1211 PetscCall(PetscFunctionListAdd(&physics, "advect", PhysicsCreate_Advect)); 1212 1213 ctx.comm = comm; 1214 ctx.cfl = 0.9; 1215 ctx.bctype = FVBC_PERIODIC; 1216 ctx.xmin = -1.0; 1217 ctx.xmax = 1.0; 1218 ctx.initial = 1; 1219 ctx.hratio = 2; 1220 maxtime = 10.0; 1221 ctx.simulation = PETSC_FALSE; 1222 PetscOptionsBegin(comm, NULL, "Finite Volume solver options", ""); 1223 PetscCall(PetscOptionsReal("-xmin", "X min", "", ctx.xmin, &ctx.xmin, NULL)); 1224 PetscCall(PetscOptionsReal("-xmax", "X max", "", ctx.xmax, &ctx.xmax, NULL)); 1225 PetscCall(PetscOptionsFList("-limit", "Name of flux imiter to use", "", limiters, lname, lname, sizeof(lname), NULL)); 1226 PetscCall(PetscOptionsFList("-physics", "Name of physics model to use", "", physics, physname, physname, sizeof(physname), NULL)); 1227 PetscCall(PetscOptionsInt("-draw", "Draw solution vector, bitwise OR of (1=initial,2=final,4=final error)", "", draw, &draw, NULL)); 1228 PetscCall(PetscOptionsString("-view_final", "Write final solution in ASCII MATLAB format to given file name", "", final_fname, final_fname, sizeof(final_fname), &view_final)); 1229 PetscCall(PetscOptionsInt("-initial", "Initial condition (depends on the physics)", "", ctx.initial, &ctx.initial, NULL)); 1230 PetscCall(PetscOptionsBool("-exact", "Compare errors with exact solution", "", ctx.exact, &ctx.exact, NULL)); 1231 PetscCall(PetscOptionsBool("-simulation", "Compare errors with reference solution", "", ctx.simulation, &ctx.simulation, NULL)); 1232 PetscCall(PetscOptionsReal("-cfl", "CFL number to time step at", "", ctx.cfl, &ctx.cfl, NULL)); 1233 PetscCall(PetscOptionsEnum("-bc_type", "Boundary condition", "", FVBCTypes, (PetscEnum)ctx.bctype, (PetscEnum *)&ctx.bctype, NULL)); 1234 PetscCall(PetscOptionsInt("-hratio", "Spacing ratio", "", ctx.hratio, &ctx.hratio, NULL)); 1235 PetscOptionsEnd(); 1236 1237 /* Choose the limiter from the list of registered limiters */ 1238 PetscCall(PetscFunctionListFind(limiters, lname, &ctx.limit2)); 1239 PetscCheck(ctx.limit2, PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "Limiter '%s' not found", lname); 1240 1241 /* Choose the physics from the list of registered models */ 1242 { 1243 PetscErrorCode (*r)(FVCtx *); 1244 PetscCall(PetscFunctionListFind(physics, physname, &r)); 1245 PetscCheck(r, PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "Physics '%s' not found", physname); 1246 /* Create the physics, will set the number of fields and their names */ 1247 PetscCall((*r)(&ctx)); 1248 } 1249 1250 /* Create a DMDA to manage the parallel grid */ 1251 PetscCall(DMDACreate1d(comm, DM_BOUNDARY_PERIODIC, 50, ctx.physics2.dof, 2, NULL, &da)); 1252 PetscCall(DMSetFromOptions(da)); 1253 PetscCall(DMSetUp(da)); 1254 /* Inform the DMDA of the field names provided by the physics. */ 1255 /* The names will be shown in the title bars when run with -ts_monitor_draw_solution */ 1256 for (i = 0; i < ctx.physics2.dof; i++) PetscCall(DMDASetFieldName(da, i, ctx.physics2.fieldname[i])); 1257 PetscCall(DMDAGetInfo(da, 0, &Mx, 0, 0, 0, 0, 0, &dof, 0, 0, 0, 0, 0)); 1258 PetscCall(DMDAGetCorners(da, &xs, 0, 0, &xm, 0, 0)); 1259 1260 /* Set coordinates of cell centers */ 1261 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)); 1262 1263 /* Allocate work space for the Finite Volume solver (so it doesn't have to be reallocated on each function evaluation) */ 1264 PetscCall(PetscMalloc4(dof * dof, &ctx.R, dof * dof, &ctx.Rinv, 2 * dof, &ctx.cjmpLR, 1 * dof, &ctx.cslope)); 1265 PetscCall(PetscMalloc3(2 * dof, &ctx.uLR, dof, &ctx.flux, dof, &ctx.speeds)); 1266 PetscCall(PetscMalloc1(2 * dof, &ctx.ub)); 1267 1268 /* Create a vector to store the solution and to save the initial state */ 1269 PetscCall(DMCreateGlobalVector(da, &X)); 1270 PetscCall(VecDuplicate(X, &X0)); 1271 PetscCall(VecDuplicate(X, &R)); 1272 1273 /* create index for slow parts and fast parts, 1274 count_slow + count_fast = Mx, counts_slow*hs = 0.5, counts_fast*hf = 0.5 */ 1275 count_slow = Mx * 3 / (3 + ctx.hratio); // compute Mx / (1.0 + ctx.hratio / 3.0); 1276 PetscCheck(count_slow % 2 == 0, PETSC_COMM_WORLD, PETSC_ERR_USER, "Please adjust grid size Mx (-da_grid_x) and hratio (-hratio) so that Mx/(1+hratio/3) is even"); 1277 count_fast = Mx - count_slow; 1278 ctx.sf = count_slow / 2; 1279 ctx.fs = ctx.sf + count_fast; 1280 PetscCall(PetscMalloc1(xm * dof, &index_slow)); 1281 PetscCall(PetscMalloc1(xm * dof, &index_fast)); 1282 PetscCall(PetscMalloc1(8 * dof, &index_slowbuffer)); 1283 ctx.lsbwidth = 4; 1284 ctx.rsbwidth = 4; 1285 1286 for (i = xs; i < xs + xm; i++) { 1287 if (i < ctx.sf - ctx.lsbwidth || i > ctx.fs + ctx.rsbwidth - 1) 1288 for (k = 0; k < dof; k++) index_slow[islow++] = i * dof + k; 1289 else if ((i >= ctx.sf - ctx.lsbwidth && i < ctx.sf) || (i > ctx.fs - 1 && i <= ctx.fs + ctx.rsbwidth - 1)) 1290 for (k = 0; k < dof; k++) index_slowbuffer[islowbuffer++] = i * dof + k; 1291 else 1292 for (k = 0; k < dof; k++) index_fast[ifast++] = i * dof + k; 1293 } 1294 PetscCall(ISCreateGeneral(PETSC_COMM_WORLD, islow, index_slow, PETSC_COPY_VALUES, &ctx.iss)); 1295 PetscCall(ISCreateGeneral(PETSC_COMM_WORLD, ifast, index_fast, PETSC_COPY_VALUES, &ctx.isf)); 1296 PetscCall(ISCreateGeneral(PETSC_COMM_WORLD, islowbuffer, index_slowbuffer, PETSC_COPY_VALUES, &ctx.issb)); 1297 1298 /* Create a time-stepping object */ 1299 PetscCall(TSCreate(comm, &ts)); 1300 PetscCall(TSSetDM(ts, da)); 1301 PetscCall(TSSetRHSFunction(ts, R, FVRHSFunction_2WaySplit, &ctx)); 1302 PetscCall(TSRHSSplitSetIS(ts, "slow", ctx.iss)); 1303 PetscCall(TSRHSSplitSetIS(ts, "slowbuffer", ctx.issb)); 1304 PetscCall(TSRHSSplitSetIS(ts, "fast", ctx.isf)); 1305 PetscCall(TSRHSSplitSetRHSFunction(ts, "slow", NULL, FVRHSFunctionslow_2WaySplit, &ctx)); 1306 PetscCall(TSRHSSplitSetRHSFunction(ts, "fast", NULL, FVRHSFunctionfast_2WaySplit, &ctx)); 1307 PetscCall(TSRHSSplitSetRHSFunction(ts, "slowbuffer", NULL, FVRHSFunctionslowbuffer_2WaySplit, &ctx)); 1308 1309 PetscCall(TSSetType(ts, TSMPRK)); 1310 PetscCall(TSSetMaxTime(ts, maxtime)); 1311 PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER)); 1312 1313 /* Compute initial conditions and starting time step */ 1314 PetscCall(FVSample_2WaySplit(&ctx, da, 0, X0)); 1315 PetscCall(FVRHSFunction_2WaySplit(ts, 0, X0, X, (void *)&ctx)); /* Initial function evaluation, only used to determine max speed */ 1316 PetscCall(VecCopy(X0, X)); /* The function value was not used so we set X=X0 again */ 1317 PetscCall(TSSetTimeStep(ts, ctx.cfl / ctx.cfl_idt)); 1318 PetscCall(TSSetFromOptions(ts)); /* Take runtime options */ 1319 PetscCall(SolutionStatsView(da, X, PETSC_VIEWER_STDOUT_WORLD)); 1320 { 1321 PetscInt steps; 1322 PetscScalar mass_initial, mass_final, mass_difference, mass_differenceg; 1323 const PetscScalar *ptr_X, *ptr_X0; 1324 const PetscReal hs = (ctx.xmax - ctx.xmin) * 3.0 / 4.0 / count_slow; 1325 const PetscReal hf = (ctx.xmax - ctx.xmin) / 4.0 / count_fast; 1326 1327 PetscCall(TSSolve(ts, X)); 1328 PetscCall(TSGetSolveTime(ts, &ptime)); 1329 PetscCall(TSGetStepNumber(ts, &steps)); 1330 /* calculate the total mass at initial time and final time */ 1331 mass_initial = 0.0; 1332 mass_final = 0.0; 1333 PetscCall(DMDAVecGetArrayRead(da, X0, (void *)&ptr_X0)); 1334 PetscCall(DMDAVecGetArrayRead(da, X, (void *)&ptr_X)); 1335 for (i = xs; i < xs + xm; i++) { 1336 if (i < ctx.sf || i > ctx.fs - 1) { 1337 for (k = 0; k < dof; k++) { 1338 mass_initial = mass_initial + hs * ptr_X0[i * dof + k]; 1339 mass_final = mass_final + hs * ptr_X[i * dof + k]; 1340 } 1341 } else { 1342 for (k = 0; k < dof; k++) { 1343 mass_initial = mass_initial + hf * ptr_X0[i * dof + k]; 1344 mass_final = mass_final + hf * ptr_X[i * dof + k]; 1345 } 1346 } 1347 } 1348 PetscCall(DMDAVecRestoreArrayRead(da, X0, (void *)&ptr_X0)); 1349 PetscCall(DMDAVecRestoreArrayRead(da, X, (void *)&ptr_X)); 1350 mass_difference = mass_final - mass_initial; 1351 PetscCall(MPIU_Allreduce(&mass_difference, &mass_differenceg, 1, MPIU_SCALAR, MPIU_SUM, comm)); 1352 PetscCall(PetscPrintf(comm, "Mass difference %g\n", (double)mass_differenceg)); 1353 PetscCall(PetscPrintf(comm, "Final time %g, steps %" PetscInt_FMT "\n", (double)ptime, steps)); 1354 PetscCall(PetscPrintf(comm, "Maximum allowable stepsize according to CFL %g\n", (double)(1.0 / ctx.cfl_idt))); 1355 if (ctx.exact) { 1356 PetscReal nrm1 = 0; 1357 PetscCall(SolutionErrorNorms_2WaySplit(&ctx, da, ptime, X, &nrm1)); 1358 PetscCall(PetscPrintf(comm, "Error ||x-x_e||_1 %g\n", (double)nrm1)); 1359 } 1360 if (ctx.simulation) { 1361 PetscReal nrm1 = 0; 1362 PetscViewer fd; 1363 char filename[PETSC_MAX_PATH_LEN] = "binaryoutput"; 1364 Vec XR; 1365 PetscBool flg; 1366 const PetscScalar *ptr_XR; 1367 PetscCall(PetscOptionsGetString(NULL, NULL, "-f", filename, sizeof(filename), &flg)); 1368 PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_USER, "Must indicate binary file with the -f option"); 1369 PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, filename, FILE_MODE_READ, &fd)); 1370 PetscCall(VecDuplicate(X0, &XR)); 1371 PetscCall(VecLoad(XR, fd)); 1372 PetscCall(PetscViewerDestroy(&fd)); 1373 PetscCall(VecGetArrayRead(X, &ptr_X)); 1374 PetscCall(VecGetArrayRead(XR, &ptr_XR)); 1375 for (i = xs; i < xs + xm; i++) { 1376 if (i < ctx.sf || i > ctx.fs - 1) 1377 for (k = 0; k < dof; k++) nrm1 = nrm1 + hs * PetscAbs(ptr_X[i * dof + k] - ptr_XR[i * dof + k]); 1378 else 1379 for (k = 0; k < dof; k++) nrm1 = nrm1 + hf * PetscAbs(ptr_X[i * dof + k] - ptr_XR[i * dof + k]); 1380 } 1381 PetscCall(VecRestoreArrayRead(X, &ptr_X)); 1382 PetscCall(VecRestoreArrayRead(XR, &ptr_XR)); 1383 PetscCall(PetscPrintf(comm, "Error ||x-x_e||_1 %g\n", (double)nrm1)); 1384 PetscCall(VecDestroy(&XR)); 1385 } 1386 } 1387 1388 PetscCall(SolutionStatsView(da, X, PETSC_VIEWER_STDOUT_WORLD)); 1389 if (draw & 0x1) PetscCall(VecView(X0, PETSC_VIEWER_DRAW_WORLD)); 1390 if (draw & 0x2) PetscCall(VecView(X, PETSC_VIEWER_DRAW_WORLD)); 1391 if (draw & 0x4) { 1392 Vec Y; 1393 PetscCall(VecDuplicate(X, &Y)); 1394 PetscCall(FVSample_2WaySplit(&ctx, da, ptime, Y)); 1395 PetscCall(VecAYPX(Y, -1, X)); 1396 PetscCall(VecView(Y, PETSC_VIEWER_DRAW_WORLD)); 1397 PetscCall(VecDestroy(&Y)); 1398 } 1399 1400 if (view_final) { 1401 PetscViewer viewer; 1402 PetscCall(PetscViewerASCIIOpen(PETSC_COMM_WORLD, final_fname, &viewer)); 1403 PetscCall(PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_MATLAB)); 1404 PetscCall(VecView(X, viewer)); 1405 PetscCall(PetscViewerPopFormat(viewer)); 1406 PetscCall(PetscViewerDestroy(&viewer)); 1407 } 1408 1409 /* Clean up */ 1410 PetscCall((*ctx.physics2.destroy)(ctx.physics2.user)); 1411 for (i = 0; i < ctx.physics2.dof; i++) PetscCall(PetscFree(ctx.physics2.fieldname[i])); 1412 PetscCall(PetscFree(ctx.physics2.bcinflowindex)); 1413 PetscCall(PetscFree(ctx.ub)); 1414 PetscCall(PetscFree4(ctx.R, ctx.Rinv, ctx.cjmpLR, ctx.cslope)); 1415 PetscCall(PetscFree3(ctx.uLR, ctx.flux, ctx.speeds)); 1416 PetscCall(VecDestroy(&X)); 1417 PetscCall(VecDestroy(&X0)); 1418 PetscCall(VecDestroy(&R)); 1419 PetscCall(DMDestroy(&da)); 1420 PetscCall(TSDestroy(&ts)); 1421 PetscCall(ISDestroy(&ctx.iss)); 1422 PetscCall(ISDestroy(&ctx.isf)); 1423 PetscCall(ISDestroy(&ctx.issb)); 1424 PetscCall(PetscFree(index_slow)); 1425 PetscCall(PetscFree(index_fast)); 1426 PetscCall(PetscFree(index_slowbuffer)); 1427 PetscCall(PetscFunctionListDestroy(&limiters)); 1428 PetscCall(PetscFunctionListDestroy(&physics)); 1429 PetscCall(PetscFinalize()); 1430 return 0; 1431 } 1432 1433 /*TEST 1434 1435 build: 1436 requires: !complex !single 1437 depends: finitevolume1d.c 1438 1439 test: 1440 suffix: 1 1441 args: -da_grid_x 60 -initial 7 -xmin -1 -xmax 1 -hratio 2 -limit mc -ts_dt 0.025 -ts_max_steps 24 -ts_type mprk -ts_mprk_type 2a22 1442 output_file: output/ex4_1.out 1443 1444 test: 1445 suffix: 2 1446 args: -da_grid_x 60 -initial 7 -xmin -1 -xmax 1 -hratio 2 -limit mc -ts_dt 0.025 -ts_max_steps 24 -ts_type mprk -ts_mprk_type 2a22 -ts_use_splitrhsfunction 1 1447 output_file: output/ex4_1.out 1448 1449 test: 1450 suffix: 3 1451 args: -da_grid_x 40 -initial 1 -hratio 2 -limit mc -ts_dt 0.1 -ts_max_steps 24 -ts_max_time 7.0 -ts_type mprk -ts_mprk_type 2a22 -physics shallow -bc_type outflow -xmin 0 -xmax 50 -ts_use_splitrhsfunction 0 1452 output_file: output/ex4_3.out 1453 1454 test: 1455 suffix: 4 1456 nsize: 2 1457 args: -da_grid_x 40 -initial 1 -hratio 2 -limit mc -ts_dt 0.1 -ts_max_steps 24 -ts_max_time 7.0 -ts_type mprk -ts_mprk_type 2a22 -physics shallow -bc_type outflow -xmin 0 -xmax 50 -ts_use_splitrhsfunction 1 1458 output_file: output/ex4_3.out 1459 1460 test: 1461 suffix: 5 1462 nsize: 4 1463 args: -da_grid_x 40 -initial 1 -hratio 2 -limit mc -ts_dt 0.1 -ts_max_steps 24 -ts_max_time 7.0 -ts_type mprk -ts_mprk_type 2a22 -physics shallow -bc_type outflow -xmin 0 -xmax 50 -ts_use_splitrhsfunction 1 1464 output_file: output/ex4_3.out 1465 TEST*/ 1466