1 /* 2 Note: 3 To performance a convergence study, a reference solution can be generated with a small stepsize and stored into a binary file, e.g. -ts_type ssp -ts_time_step 1e-5 -ts_max_steps 30000 -ts_view_solution binary:reference.bin 4 Errors can be computed in the following runs with -simulation -f reference.bin 5 6 Multirate RK options: 7 -ts_rk_dtratio is the ratio between larger time step size and small time step size 8 -ts_rk_multirate_type has three choices: none (for single step RK) 9 combined (for multirate method and user just need to provide the same RHS with the single step RK) 10 partitioned (for multiraet method and user need to provide two RHS, one is for fast components and the orther is for slow components 11 */ 12 13 static const char help[] = "1D periodic Finite Volume solver in slope-limiter form with semidiscrete time stepping.\n" 14 " advection - Variable coefficient scalar advection\n" 15 " u_t + (a*u)_x = 0\n" 16 " a is a piecewise function with two values say a0 and a1 (both a0 and a1 are constant).\n" 17 " in this toy problem we assume a=a0 when x<0 and a=a1 when x>0 with a0<a1 which has the same discontinuous point with initial condition.\n" 18 " we don't provide an exact solution, so you need to generate reference solution to do convergnece staudy,\n" 19 " more precisely, you need firstly generate a reference to a binary file say file.bin, then on command line,\n" 20 " you should type -simulation -f file.bin.\n" 21 " you can choose the number of grids by -da_grid_x.\n" 22 " you can choose the value of a by -physics_advect_a1 and -physics_advect_a2.\n"; 23 24 #include <petscts.h> 25 #include <petscdm.h> 26 #include <petscdmda.h> 27 #include <petscdraw.h> 28 #include <petsc/private/tsimpl.h> 29 30 #include <petsc/private/kernels/blockinvert.h> /* For the Kernel_*_gets_* stuff for BAIJ */ 31 32 #include "finitevolume1d.h" 33 34 static inline PetscReal RangeMod(PetscReal a, PetscReal xmin, PetscReal xmax) 35 { 36 PetscReal range = xmax - xmin; 37 return xmin + PetscFmodReal(range + PetscFmodReal(a, range), range); 38 } 39 40 /* --------------------------------- Advection ----------------------------------- */ 41 42 typedef struct { 43 PetscReal a[2]; /* advective velocity */ 44 } AdvectCtx; 45 46 static PetscErrorCode PhysicsRiemann_Advect(void *vctx, PetscInt m, const PetscScalar *uL, const PetscScalar *uR, PetscScalar *flux, PetscReal *maxspeed, PetscReal x, PetscReal xmin, PetscReal xmax) 47 { 48 AdvectCtx *ctx = (AdvectCtx *)vctx; 49 PetscReal *speed; 50 51 PetscFunctionBeginUser; 52 speed = ctx->a; 53 if (x == 0 || x == xmin || x == xmax) 54 flux[0] = PetscMax(0, (speed[0] + speed[1]) / 2.0) * uL[0] + PetscMin(0, (speed[0] + speed[1]) / 2.0) * uR[0]; /* if condition need to be changed base on different problem, '0' is the discontinuous point of a */ 55 else if (x < 0) flux[0] = PetscMax(0, speed[0]) * uL[0] + PetscMin(0, speed[0]) * uR[0]; /* else if condition need to be changed based on different problem, 'x = 0' is discontinuous point of a */ 56 else flux[0] = PetscMax(0, speed[1]) * uL[0] + PetscMin(0, speed[1]) * uR[0]; 57 *maxspeed = *speed; 58 PetscFunctionReturn(PETSC_SUCCESS); 59 } 60 61 static PetscErrorCode PhysicsCharacteristic_Advect(void *vctx, PetscInt m, const PetscScalar *u, PetscScalar *X, PetscScalar *Xi, PetscReal *speeds, PetscReal x) 62 { 63 AdvectCtx *ctx = (AdvectCtx *)vctx; 64 65 PetscFunctionBeginUser; 66 X[0] = 1.; 67 Xi[0] = 1.; 68 if (x < 0) speeds[0] = ctx->a[0]; /* x is discontinuous point of a */ 69 else speeds[0] = ctx->a[1]; 70 PetscFunctionReturn(PETSC_SUCCESS); 71 } 72 73 static PetscErrorCode PhysicsSample_Advect(void *vctx, PetscInt initial, FVBCType bctype, PetscReal xmin, PetscReal xmax, PetscReal t, PetscReal x, PetscReal *u) 74 { 75 AdvectCtx *ctx = (AdvectCtx *)vctx; 76 PetscReal *a = ctx->a, x0; 77 78 PetscFunctionBeginUser; 79 if (x < 0) { /* x is cell center */ 80 switch (bctype) { 81 case FVBC_OUTFLOW: 82 x0 = x - a[0] * t; 83 break; 84 case FVBC_PERIODIC: 85 x0 = RangeMod(x - a[0] * t, xmin, xmax); 86 break; 87 default: 88 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "unknown BCType"); 89 } 90 switch (initial) { 91 case 0: 92 u[0] = (x0 < 0) ? 1 : -1; 93 break; 94 case 1: 95 u[0] = (x0 < 0) ? -1 : 1; 96 break; 97 case 2: 98 u[0] = (0 < x0 && x0 < 1) ? 1 : 0; 99 break; 100 case 3: 101 u[0] = PetscSinReal(2 * PETSC_PI * x0); 102 break; 103 case 4: 104 u[0] = PetscAbs(x0); 105 break; 106 case 5: 107 u[0] = (x0 < 0 || x0 > 0.5) ? 0 : PetscSqr(PetscSinReal(2 * PETSC_PI * x0)); 108 break; 109 case 6: 110 u[0] = (x0 < 0) ? 0 : ((x0 < 1) ? x0 : ((x0 < 2) ? 2 - x0 : 0)); 111 break; 112 case 7: 113 u[0] = PetscPowReal(PetscSinReal(PETSC_PI * x0), 10.0); 114 break; 115 default: 116 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "unknown initial condition"); 117 } 118 } else { 119 switch (bctype) { 120 case FVBC_OUTFLOW: 121 x0 = x - a[1] * t; 122 break; 123 case FVBC_PERIODIC: 124 x0 = RangeMod(x - a[1] * t, xmin, xmax); 125 break; 126 default: 127 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "unknown BCType"); 128 } 129 switch (initial) { 130 case 0: 131 u[0] = (x0 < 0) ? 1 : -1; 132 break; 133 case 1: 134 u[0] = (x0 < 0) ? -1 : 1; 135 break; 136 case 2: 137 u[0] = (0 < x0 && x0 < 1) ? 1 : 0; 138 break; 139 case 3: 140 u[0] = PetscSinReal(2 * PETSC_PI * x0); 141 break; 142 case 4: 143 u[0] = PetscAbs(x0); 144 break; 145 case 5: 146 u[0] = (x0 < 0 || x0 > 0.5) ? 0 : PetscSqr(PetscSinReal(2 * PETSC_PI * x0)); 147 break; 148 case 6: 149 u[0] = (x0 < 0) ? 0 : ((x0 < 1) ? x0 : ((x0 < 2) ? 2 - x0 : 0)); 150 break; 151 case 7: 152 u[0] = PetscPowReal(PetscSinReal(PETSC_PI * x0), 10.0); 153 break; 154 default: 155 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "unknown initial condition"); 156 } 157 } 158 PetscFunctionReturn(PETSC_SUCCESS); 159 } 160 161 static PetscErrorCode PhysicsCreate_Advect(FVCtx *ctx) 162 { 163 AdvectCtx *user; 164 165 PetscFunctionBeginUser; 166 PetscCall(PetscNew(&user)); 167 ctx->physics.sample = PhysicsSample_Advect; 168 ctx->physics.riemann = PhysicsRiemann_Advect; 169 ctx->physics.characteristic = PhysicsCharacteristic_Advect; 170 ctx->physics.destroy = PhysicsDestroy_SimpleFree; 171 ctx->physics.user = user; 172 ctx->physics.dof = 1; 173 PetscCall(PetscStrallocpy("u", &ctx->physics.fieldname[0])); 174 user->a[0] = 1; 175 user->a[1] = 1; 176 PetscOptionsBegin(ctx->comm, ctx->prefix, "Options for advection", ""); 177 { 178 PetscCall(PetscOptionsReal("-physics_advect_a1", "Speed1", "", user->a[0], &user->a[0], NULL)); 179 PetscCall(PetscOptionsReal("-physics_advect_a2", "Speed2", "", user->a[1], &user->a[1], NULL)); 180 } 181 PetscOptionsEnd(); 182 PetscFunctionReturn(PETSC_SUCCESS); 183 } 184 185 /* --------------------------------- Finite Volume Solver for slow parts ----------------------------------- */ 186 187 PetscErrorCode FVRHSFunctionslow(TS ts, PetscReal time, Vec X, Vec F, void *vctx) 188 { 189 FVCtx *ctx = (FVCtx *)vctx; 190 PetscInt i, j, k, Mx, dof, xs, xm, len_slow; 191 PetscReal hx, cfl_idt = 0; 192 PetscScalar *x, *f, *slope; 193 Vec Xloc; 194 DM da; 195 196 PetscFunctionBeginUser; 197 PetscCall(TSGetDM(ts, &da)); 198 PetscCall(DMGetLocalVector(da, &Xloc)); 199 PetscCall(DMDAGetInfo(da, 0, &Mx, 0, 0, 0, 0, 0, &dof, 0, 0, 0, 0, 0)); 200 hx = (ctx->xmax - ctx->xmin) / Mx; 201 PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, Xloc)); 202 PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, Xloc)); 203 PetscCall(VecZeroEntries(F)); 204 PetscCall(DMDAVecGetArray(da, Xloc, &x)); 205 PetscCall(VecGetArray(F, &f)); 206 PetscCall(DMDAGetArray(da, PETSC_TRUE, &slope)); 207 PetscCall(DMDAGetCorners(da, &xs, 0, 0, &xm, 0, 0)); 208 PetscCall(ISGetSize(ctx->iss, &len_slow)); 209 210 if (ctx->bctype == FVBC_OUTFLOW) { 211 for (i = xs - 2; i < 0; i++) { 212 for (j = 0; j < dof; j++) x[i * dof + j] = x[j]; 213 } 214 for (i = Mx; i < xs + xm + 2; i++) { 215 for (j = 0; j < dof; j++) x[i * dof + j] = x[(xs + xm - 1) * dof + j]; 216 } 217 } 218 for (i = xs - 1; i < xs + xm + 1; i++) { 219 struct _LimitInfo info; 220 PetscScalar *cjmpL, *cjmpR; 221 if (i < len_slow + 1) { 222 /* Determine the right eigenvectors R, where A = R \Lambda R^{-1} */ 223 PetscCall((*ctx->physics.characteristic)(ctx->physics.user, dof, &x[i * dof], ctx->R, ctx->Rinv, ctx->speeds, ctx->xmin + hx * i)); 224 /* Evaluate jumps across interfaces (i-1, i) and (i, i+1), put in characteristic basis */ 225 PetscCall(PetscArrayzero(ctx->cjmpLR, 2 * dof)); 226 cjmpL = &ctx->cjmpLR[0]; 227 cjmpR = &ctx->cjmpLR[dof]; 228 for (j = 0; j < dof; j++) { 229 PetscScalar jmpL, jmpR; 230 jmpL = x[(i + 0) * dof + j] - x[(i - 1) * dof + j]; 231 jmpR = x[(i + 1) * dof + j] - x[(i + 0) * dof + j]; 232 for (k = 0; k < dof; k++) { 233 cjmpL[k] += ctx->Rinv[k + j * dof] * jmpL; 234 cjmpR[k] += ctx->Rinv[k + j * dof] * jmpR; 235 } 236 } 237 /* Apply limiter to the left and right characteristic jumps */ 238 info.m = dof; 239 info.hx = hx; 240 (*ctx->limit)(&info, cjmpL, cjmpR, ctx->cslope); 241 for (j = 0; j < dof; j++) ctx->cslope[j] /= hx; /* rescale to a slope */ 242 for (j = 0; j < dof; j++) { 243 PetscScalar tmp = 0; 244 for (k = 0; k < dof; k++) tmp += ctx->R[j + k * dof] * ctx->cslope[k]; 245 slope[i * dof + j] = tmp; 246 } 247 } 248 } 249 250 for (i = xs; i < xs + xm + 1; i++) { 251 PetscReal maxspeed; 252 PetscScalar *uL, *uR; 253 if (i < len_slow) { /* slow parts can be changed based on a */ 254 uL = &ctx->uLR[0]; 255 uR = &ctx->uLR[dof]; 256 for (j = 0; j < dof; j++) { 257 uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hx / 2; 258 uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hx / 2; 259 } 260 PetscCall((*ctx->physics.riemann)(ctx->physics.user, dof, uL, uR, ctx->flux, &maxspeed, ctx->xmin + hx * i, ctx->xmin, ctx->xmax)); 261 cfl_idt = PetscMax(cfl_idt, PetscAbsScalar(maxspeed / hx)); /* Max allowable value of 1/Delta t */ 262 if (i > xs) { 263 for (j = 0; j < dof; j++) f[(i - 1) * dof + j] -= ctx->flux[j] / hx; 264 } 265 if (i < xs + xm) { 266 for (j = 0; j < dof; j++) f[i * dof + j] += ctx->flux[j] / hx; 267 } 268 } 269 if (i == len_slow) { /* slow parts can be changed based on a */ 270 uL = &ctx->uLR[0]; 271 uR = &ctx->uLR[dof]; 272 for (j = 0; j < dof; j++) { 273 uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hx / 2; 274 uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hx / 2; 275 } 276 PetscCall((*ctx->physics.riemann)(ctx->physics.user, dof, uL, uR, ctx->flux, &maxspeed, ctx->xmin + hx * i, ctx->xmin, ctx->xmax)); 277 cfl_idt = PetscMax(cfl_idt, PetscAbsScalar(maxspeed / hx)); /* Max allowable value of 1/Delta t */ 278 if (i > xs) { 279 for (j = 0; j < dof; j++) f[(i - 1) * dof + j] -= ctx->flux[j] / hx; 280 } 281 } 282 } 283 PetscCall(DMDAVecRestoreArray(da, Xloc, &x)); 284 PetscCall(VecRestoreArray(F, &f)); 285 PetscCall(DMDARestoreArray(da, PETSC_TRUE, &slope)); 286 PetscCall(DMRestoreLocalVector(da, &Xloc)); 287 PetscFunctionReturn(PETSC_SUCCESS); 288 } 289 290 /* --------------------------------- Finite Volume Solver for fast parts ----------------------------------- */ 291 292 PetscErrorCode FVRHSFunctionfast(TS ts, PetscReal time, Vec X, Vec F, void *vctx) 293 { 294 FVCtx *ctx = (FVCtx *)vctx; 295 PetscInt i, j, k, Mx, dof, xs, xm, len_slow; 296 PetscReal hx, cfl_idt = 0; 297 PetscScalar *x, *f, *slope; 298 Vec Xloc; 299 DM da; 300 301 PetscFunctionBeginUser; 302 PetscCall(TSGetDM(ts, &da)); 303 PetscCall(DMGetLocalVector(da, &Xloc)); 304 PetscCall(DMDAGetInfo(da, 0, &Mx, 0, 0, 0, 0, 0, &dof, 0, 0, 0, 0, 0)); 305 hx = (ctx->xmax - ctx->xmin) / Mx; 306 PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, Xloc)); 307 PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, Xloc)); 308 PetscCall(VecZeroEntries(F)); 309 PetscCall(DMDAVecGetArray(da, Xloc, &x)); 310 PetscCall(VecGetArray(F, &f)); 311 PetscCall(DMDAGetArray(da, PETSC_TRUE, &slope)); 312 PetscCall(DMDAGetCorners(da, &xs, 0, 0, &xm, 0, 0)); 313 PetscCall(ISGetSize(ctx->iss, &len_slow)); 314 315 if (ctx->bctype == FVBC_OUTFLOW) { 316 for (i = xs - 2; i < 0; i++) { 317 for (j = 0; j < dof; j++) x[i * dof + j] = x[j]; 318 } 319 for (i = Mx; i < xs + xm + 2; i++) { 320 for (j = 0; j < dof; j++) x[i * dof + j] = x[(xs + xm - 1) * dof + j]; 321 } 322 } 323 for (i = xs - 1; i < xs + xm + 1; i++) { 324 struct _LimitInfo info; 325 PetscScalar *cjmpL, *cjmpR; 326 if (i > len_slow - 2) { 327 PetscCall((*ctx->physics.characteristic)(ctx->physics.user, dof, &x[i * dof], ctx->R, ctx->Rinv, ctx->speeds, ctx->xmin + hx * i)); 328 PetscCall(PetscArrayzero(ctx->cjmpLR, 2 * dof)); 329 cjmpL = &ctx->cjmpLR[0]; 330 cjmpR = &ctx->cjmpLR[dof]; 331 for (j = 0; j < dof; j++) { 332 PetscScalar jmpL, jmpR; 333 jmpL = x[(i + 0) * dof + j] - x[(i - 1) * dof + j]; 334 jmpR = x[(i + 1) * dof + j] - x[(i + 0) * dof + j]; 335 for (k = 0; k < dof; k++) { 336 cjmpL[k] += ctx->Rinv[k + j * dof] * jmpL; 337 cjmpR[k] += ctx->Rinv[k + j * dof] * jmpR; 338 } 339 } 340 /* Apply limiter to the left and right characteristic jumps */ 341 info.m = dof; 342 info.hx = hx; 343 (*ctx->limit)(&info, cjmpL, cjmpR, ctx->cslope); 344 for (j = 0; j < dof; j++) ctx->cslope[j] /= hx; /* rescale to a slope */ 345 for (j = 0; j < dof; j++) { 346 PetscScalar tmp = 0; 347 for (k = 0; k < dof; k++) tmp += ctx->R[j + k * dof] * ctx->cslope[k]; 348 slope[i * dof + j] = tmp; 349 } 350 } 351 } 352 353 for (i = xs; i < xs + xm + 1; i++) { 354 PetscReal maxspeed; 355 PetscScalar *uL, *uR; 356 if (i > len_slow) { 357 uL = &ctx->uLR[0]; 358 uR = &ctx->uLR[dof]; 359 for (j = 0; j < dof; j++) { 360 uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hx / 2; 361 uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hx / 2; 362 } 363 PetscCall((*ctx->physics.riemann)(ctx->physics.user, dof, uL, uR, ctx->flux, &maxspeed, ctx->xmin + hx * i, ctx->xmin, ctx->xmax)); 364 cfl_idt = PetscMax(cfl_idt, PetscAbsScalar(maxspeed / hx)); /* Max allowable value of 1/Delta t */ 365 if (i > xs) { 366 for (j = 0; j < dof; j++) f[(i - len_slow - 1) * dof + j] -= ctx->flux[j] / hx; 367 } 368 if (i < xs + xm) { 369 for (j = 0; j < dof; j++) f[(i - len_slow) * dof + j] += ctx->flux[j] / hx; 370 } 371 } 372 if (i == len_slow) { 373 uL = &ctx->uLR[0]; 374 uR = &ctx->uLR[dof]; 375 for (j = 0; j < dof; j++) { 376 uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hx / 2; 377 uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hx / 2; 378 } 379 PetscCall((*ctx->physics.riemann)(ctx->physics.user, dof, uL, uR, ctx->flux, &maxspeed, ctx->xmin + hx * i, ctx->xmin, ctx->xmax)); 380 cfl_idt = PetscMax(cfl_idt, PetscAbsScalar(maxspeed / hx)); /* Max allowable value of 1/Delta t */ 381 if (i < xs + xm) { 382 for (j = 0; j < dof; j++) f[(i - len_slow) * dof + j] += ctx->flux[j] / hx; 383 } 384 } 385 } 386 PetscCall(DMDAVecRestoreArray(da, Xloc, &x)); 387 PetscCall(VecRestoreArray(F, &f)); 388 PetscCall(DMDARestoreArray(da, PETSC_TRUE, &slope)); 389 PetscCall(DMRestoreLocalVector(da, &Xloc)); 390 PetscFunctionReturn(PETSC_SUCCESS); 391 } 392 393 PetscErrorCode FVRHSFunctionslow2(TS ts, PetscReal time, Vec X, Vec F, void *vctx) 394 { 395 FVCtx *ctx = (FVCtx *)vctx; 396 PetscInt i, j, k, Mx, dof, xs, xm, len_slow1, len_slow2; 397 PetscReal hx, cfl_idt = 0; 398 PetscScalar *x, *f, *slope; 399 Vec Xloc; 400 DM da; 401 402 PetscFunctionBeginUser; 403 PetscCall(TSGetDM(ts, &da)); 404 PetscCall(DMGetLocalVector(da, &Xloc)); 405 PetscCall(DMDAGetInfo(da, 0, &Mx, 0, 0, 0, 0, 0, &dof, 0, 0, 0, 0, 0)); 406 hx = (ctx->xmax - ctx->xmin) / Mx; 407 PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, Xloc)); 408 PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, Xloc)); 409 PetscCall(VecZeroEntries(F)); 410 PetscCall(DMDAVecGetArray(da, Xloc, &x)); 411 PetscCall(VecGetArray(F, &f)); 412 PetscCall(DMDAGetArray(da, PETSC_TRUE, &slope)); 413 PetscCall(DMDAGetCorners(da, &xs, 0, 0, &xm, 0, 0)); 414 PetscCall(ISGetSize(ctx->iss, &len_slow1)); 415 PetscCall(ISGetSize(ctx->iss2, &len_slow2)); 416 if (ctx->bctype == FVBC_OUTFLOW) { 417 for (i = xs - 2; i < 0; i++) { 418 for (j = 0; j < dof; j++) x[i * dof + j] = x[j]; 419 } 420 for (i = Mx; i < xs + xm + 2; i++) { 421 for (j = 0; j < dof; j++) x[i * dof + j] = x[(xs + xm - 1) * dof + j]; 422 } 423 } 424 for (i = xs - 1; i < xs + xm + 1; i++) { 425 struct _LimitInfo info; 426 PetscScalar *cjmpL, *cjmpR; 427 if (i < len_slow1 + len_slow2 + 1 && i > len_slow1 - 2) { 428 PetscCall((*ctx->physics.characteristic)(ctx->physics.user, dof, &x[i * dof], ctx->R, ctx->Rinv, ctx->speeds, ctx->xmin + hx * i)); 429 PetscCall(PetscArrayzero(ctx->cjmpLR, 2 * dof)); 430 cjmpL = &ctx->cjmpLR[0]; 431 cjmpR = &ctx->cjmpLR[dof]; 432 for (j = 0; j < dof; j++) { 433 PetscScalar jmpL, jmpR; 434 jmpL = x[(i + 0) * dof + j] - x[(i - 1) * dof + j]; 435 jmpR = x[(i + 1) * dof + j] - x[(i + 0) * dof + j]; 436 for (k = 0; k < dof; k++) { 437 cjmpL[k] += ctx->Rinv[k + j * dof] * jmpL; 438 cjmpR[k] += ctx->Rinv[k + j * dof] * jmpR; 439 } 440 } 441 /* Apply limiter to the left and right characteristic jumps */ 442 info.m = dof; 443 info.hx = hx; 444 (*ctx->limit)(&info, cjmpL, cjmpR, ctx->cslope); 445 for (j = 0; j < dof; j++) ctx->cslope[j] /= hx; /* rescale to a slope */ 446 for (j = 0; j < dof; j++) { 447 PetscScalar tmp = 0; 448 for (k = 0; k < dof; k++) tmp += ctx->R[j + k * dof] * ctx->cslope[k]; 449 slope[i * dof + j] = tmp; 450 } 451 } 452 } 453 454 for (i = xs; i < xs + xm + 1; i++) { 455 PetscReal maxspeed; 456 PetscScalar *uL, *uR; 457 if (i < len_slow1 + len_slow2 && i > len_slow1) { /* slow parts can be changed based on a */ 458 uL = &ctx->uLR[0]; 459 uR = &ctx->uLR[dof]; 460 for (j = 0; j < dof; j++) { 461 uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hx / 2; 462 uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hx / 2; 463 } 464 PetscCall((*ctx->physics.riemann)(ctx->physics.user, dof, uL, uR, ctx->flux, &maxspeed, ctx->xmin + hx * i, ctx->xmin, ctx->xmax)); 465 cfl_idt = PetscMax(cfl_idt, PetscAbsScalar(maxspeed / hx)); /* Max allowable value of 1/Delta t */ 466 if (i > xs) { 467 for (j = 0; j < dof; j++) f[(i - len_slow1 - 1) * dof + j] -= ctx->flux[j] / hx; 468 } 469 if (i < xs + xm) { 470 for (j = 0; j < dof; j++) f[(i - len_slow1) * dof + j] += ctx->flux[j] / hx; 471 } 472 } 473 if (i == len_slow1 + len_slow2) { /* slow parts can be changed based on a */ 474 uL = &ctx->uLR[0]; 475 uR = &ctx->uLR[dof]; 476 for (j = 0; j < dof; j++) { 477 uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hx / 2; 478 uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hx / 2; 479 } 480 PetscCall((*ctx->physics.riemann)(ctx->physics.user, dof, uL, uR, ctx->flux, &maxspeed, ctx->xmin + hx * i, ctx->xmin, ctx->xmax)); 481 cfl_idt = PetscMax(cfl_idt, PetscAbsScalar(maxspeed / hx)); /* Max allowable value of 1/Delta t */ 482 if (i > xs) { 483 for (j = 0; j < dof; j++) f[(i - len_slow1 - 1) * dof + j] -= ctx->flux[j] / hx; 484 } 485 } 486 if (i == len_slow1) { /* slow parts can be changed based on a */ 487 uL = &ctx->uLR[0]; 488 uR = &ctx->uLR[dof]; 489 for (j = 0; j < dof; j++) { 490 uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hx / 2; 491 uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hx / 2; 492 } 493 PetscCall((*ctx->physics.riemann)(ctx->physics.user, dof, uL, uR, ctx->flux, &maxspeed, ctx->xmin + hx * i, ctx->xmin, ctx->xmax)); 494 cfl_idt = PetscMax(cfl_idt, PetscAbsScalar(maxspeed / hx)); /* Max allowable value of 1/Delta t */ 495 if (i < xs + xm) { 496 for (j = 0; j < dof; j++) f[(i - len_slow1) * dof + j] += ctx->flux[j] / hx; 497 } 498 } 499 } 500 501 PetscCall(DMDAVecRestoreArray(da, Xloc, &x)); 502 PetscCall(VecRestoreArray(F, &f)); 503 PetscCall(DMDARestoreArray(da, PETSC_TRUE, &slope)); 504 PetscCall(DMRestoreLocalVector(da, &Xloc)); 505 PetscFunctionReturn(PETSC_SUCCESS); 506 } 507 508 PetscErrorCode FVRHSFunctionfast2(TS ts, PetscReal time, Vec X, Vec F, void *vctx) 509 { 510 FVCtx *ctx = (FVCtx *)vctx; 511 PetscInt i, j, k, Mx, dof, xs, xm, len_slow1, len_slow2; 512 PetscReal hx, cfl_idt = 0; 513 PetscScalar *x, *f, *slope; 514 Vec Xloc; 515 DM da; 516 517 PetscFunctionBeginUser; 518 PetscCall(TSGetDM(ts, &da)); 519 PetscCall(DMGetLocalVector(da, &Xloc)); 520 PetscCall(DMDAGetInfo(da, 0, &Mx, 0, 0, 0, 0, 0, &dof, 0, 0, 0, 0, 0)); 521 hx = (ctx->xmax - ctx->xmin) / Mx; 522 PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, Xloc)); 523 PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, Xloc)); 524 PetscCall(VecZeroEntries(F)); 525 PetscCall(DMDAVecGetArray(da, Xloc, &x)); 526 PetscCall(VecGetArray(F, &f)); 527 PetscCall(DMDAGetArray(da, PETSC_TRUE, &slope)); 528 PetscCall(DMDAGetCorners(da, &xs, 0, 0, &xm, 0, 0)); 529 PetscCall(ISGetSize(ctx->iss, &len_slow1)); 530 PetscCall(ISGetSize(ctx->iss2, &len_slow2)); 531 532 if (ctx->bctype == FVBC_OUTFLOW) { 533 for (i = xs - 2; i < 0; i++) { 534 for (j = 0; j < dof; j++) x[i * dof + j] = x[j]; 535 } 536 for (i = Mx; i < xs + xm + 2; i++) { 537 for (j = 0; j < dof; j++) x[i * dof + j] = x[(xs + xm - 1) * dof + j]; 538 } 539 } 540 for (i = xs - 1; i < xs + xm + 1; i++) { 541 struct _LimitInfo info; 542 PetscScalar *cjmpL, *cjmpR; 543 if (i > len_slow1 + len_slow2 - 2) { 544 PetscCall((*ctx->physics.characteristic)(ctx->physics.user, dof, &x[i * dof], ctx->R, ctx->Rinv, ctx->speeds, ctx->xmin + hx * i)); 545 PetscCall(PetscArrayzero(ctx->cjmpLR, 2 * dof)); 546 cjmpL = &ctx->cjmpLR[0]; 547 cjmpR = &ctx->cjmpLR[dof]; 548 for (j = 0; j < dof; j++) { 549 PetscScalar jmpL, jmpR; 550 jmpL = x[(i + 0) * dof + j] - x[(i - 1) * dof + j]; 551 jmpR = x[(i + 1) * dof + j] - x[(i + 0) * dof + j]; 552 for (k = 0; k < dof; k++) { 553 cjmpL[k] += ctx->Rinv[k + j * dof] * jmpL; 554 cjmpR[k] += ctx->Rinv[k + j * dof] * jmpR; 555 } 556 } 557 /* Apply limiter to the left and right characteristic jumps */ 558 info.m = dof; 559 info.hx = hx; 560 (*ctx->limit)(&info, cjmpL, cjmpR, ctx->cslope); 561 for (j = 0; j < dof; j++) ctx->cslope[j] /= hx; /* rescale to a slope */ 562 for (j = 0; j < dof; j++) { 563 PetscScalar tmp = 0; 564 for (k = 0; k < dof; k++) tmp += ctx->R[j + k * dof] * ctx->cslope[k]; 565 slope[i * dof + j] = tmp; 566 } 567 } 568 } 569 570 for (i = xs; i < xs + xm + 1; i++) { 571 PetscReal maxspeed; 572 PetscScalar *uL, *uR; 573 if (i > len_slow1 + len_slow2) { 574 uL = &ctx->uLR[0]; 575 uR = &ctx->uLR[dof]; 576 for (j = 0; j < dof; j++) { 577 uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hx / 2; 578 uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hx / 2; 579 } 580 PetscCall((*ctx->physics.riemann)(ctx->physics.user, dof, uL, uR, ctx->flux, &maxspeed, ctx->xmin + hx * i, ctx->xmin, ctx->xmax)); 581 cfl_idt = PetscMax(cfl_idt, PetscAbsScalar(maxspeed / hx)); /* Max allowable value of 1/Delta t */ 582 if (i > xs) { 583 for (j = 0; j < dof; j++) f[(i - len_slow1 - len_slow2 - 1) * dof + j] -= ctx->flux[j] / hx; 584 } 585 if (i < xs + xm) { 586 for (j = 0; j < dof; j++) f[(i - len_slow1 - len_slow2) * dof + j] += ctx->flux[j] / hx; 587 } 588 } 589 if (i == len_slow1 + len_slow2) { 590 uL = &ctx->uLR[0]; 591 uR = &ctx->uLR[dof]; 592 for (j = 0; j < dof; j++) { 593 uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hx / 2; 594 uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hx / 2; 595 } 596 PetscCall((*ctx->physics.riemann)(ctx->physics.user, dof, uL, uR, ctx->flux, &maxspeed, ctx->xmin + hx * i, ctx->xmin, ctx->xmax)); 597 cfl_idt = PetscMax(cfl_idt, PetscAbsScalar(maxspeed / hx)); /* Max allowable value of 1/Delta t */ 598 if (i < xs + xm) { 599 for (j = 0; j < dof; j++) f[(i - len_slow1 - len_slow2) * dof + j] += ctx->flux[j] / hx; 600 } 601 } 602 } 603 PetscCall(DMDAVecRestoreArray(da, Xloc, &x)); 604 PetscCall(VecRestoreArray(F, &f)); 605 PetscCall(DMDARestoreArray(da, PETSC_TRUE, &slope)); 606 PetscCall(DMRestoreLocalVector(da, &Xloc)); 607 PetscFunctionReturn(PETSC_SUCCESS); 608 } 609 610 int main(int argc, char *argv[]) 611 { 612 char lname[256] = "mc", physname[256] = "advect", final_fname[256] = "solution.m"; 613 PetscFunctionList limiters = 0, physics = 0; 614 MPI_Comm comm; 615 TS ts; 616 DM da; 617 Vec X, X0, R; 618 FVCtx ctx; 619 PetscInt i, k, dof, xs, xm, Mx, draw = 0, *index_slow, *index_fast, islow = 0, ifast = 0; 620 PetscBool view_final = PETSC_FALSE; 621 PetscReal ptime; 622 623 PetscFunctionBeginUser; 624 PetscCall(PetscInitialize(&argc, &argv, 0, help)); 625 comm = PETSC_COMM_WORLD; 626 PetscCall(PetscMemzero(&ctx, sizeof(ctx))); 627 628 /* Register limiters to be available on the command line */ 629 PetscCall(PetscFunctionListAdd(&limiters, "upwind", Limit_Upwind)); 630 PetscCall(PetscFunctionListAdd(&limiters, "lax-wendroff", Limit_LaxWendroff)); 631 PetscCall(PetscFunctionListAdd(&limiters, "beam-warming", Limit_BeamWarming)); 632 PetscCall(PetscFunctionListAdd(&limiters, "fromm", Limit_Fromm)); 633 PetscCall(PetscFunctionListAdd(&limiters, "minmod", Limit_Minmod)); 634 PetscCall(PetscFunctionListAdd(&limiters, "superbee", Limit_Superbee)); 635 PetscCall(PetscFunctionListAdd(&limiters, "mc", Limit_MC)); 636 PetscCall(PetscFunctionListAdd(&limiters, "vanleer", Limit_VanLeer)); 637 PetscCall(PetscFunctionListAdd(&limiters, "vanalbada", Limit_VanAlbada)); 638 PetscCall(PetscFunctionListAdd(&limiters, "vanalbadatvd", Limit_VanAlbadaTVD)); 639 PetscCall(PetscFunctionListAdd(&limiters, "koren", Limit_Koren)); 640 PetscCall(PetscFunctionListAdd(&limiters, "korensym", Limit_KorenSym)); 641 PetscCall(PetscFunctionListAdd(&limiters, "koren3", Limit_Koren3)); 642 PetscCall(PetscFunctionListAdd(&limiters, "cada-torrilhon2", Limit_CadaTorrilhon2)); 643 PetscCall(PetscFunctionListAdd(&limiters, "cada-torrilhon3-r0p1", Limit_CadaTorrilhon3R0p1)); 644 PetscCall(PetscFunctionListAdd(&limiters, "cada-torrilhon3-r1", Limit_CadaTorrilhon3R1)); 645 PetscCall(PetscFunctionListAdd(&limiters, "cada-torrilhon3-r10", Limit_CadaTorrilhon3R10)); 646 PetscCall(PetscFunctionListAdd(&limiters, "cada-torrilhon3-r100", Limit_CadaTorrilhon3R100)); 647 648 /* Register physical models to be available on the command line */ 649 PetscCall(PetscFunctionListAdd(&physics, "advect", PhysicsCreate_Advect)); 650 651 ctx.comm = comm; 652 ctx.cfl = 0.9; 653 ctx.bctype = FVBC_PERIODIC; 654 ctx.xmin = -1.0; 655 ctx.xmax = 1.0; 656 PetscOptionsBegin(comm, NULL, "Finite Volume solver options", ""); 657 PetscCall(PetscOptionsReal("-xmin", "X min", "", ctx.xmin, &ctx.xmin, NULL)); 658 PetscCall(PetscOptionsReal("-xmax", "X max", "", ctx.xmax, &ctx.xmax, NULL)); 659 PetscCall(PetscOptionsFList("-limit", "Name of flux limiter to use", "", limiters, lname, lname, sizeof(lname), NULL)); 660 PetscCall(PetscOptionsInt("-draw", "Draw solution vector, bitwise OR of (1=initial,2=final,4=final error)", "", draw, &draw, NULL)); 661 PetscCall(PetscOptionsString("-view_final", "Write final solution in ASCII MATLAB format to given file name", "", final_fname, final_fname, sizeof(final_fname), &view_final)); 662 PetscCall(PetscOptionsInt("-initial", "Initial condition (depends on the physics)", "", ctx.initial, &ctx.initial, NULL)); 663 PetscCall(PetscOptionsBool("-simulation", "Compare errors with reference solution", "", ctx.simulation, &ctx.simulation, NULL)); 664 PetscCall(PetscOptionsReal("-cfl", "CFL number to time step at", "", ctx.cfl, &ctx.cfl, NULL)); 665 PetscCall(PetscOptionsEnum("-bc_type", "Boundary condition", "", FVBCTypes, (PetscEnum)ctx.bctype, (PetscEnum *)&ctx.bctype, NULL)); 666 PetscCall(PetscOptionsBool("-recursive_split", "Split the domain recursively", "", ctx.recursive, &ctx.recursive, NULL)); 667 PetscOptionsEnd(); 668 669 /* Choose the limiter from the list of registered limiters */ 670 PetscCall(PetscFunctionListFind(limiters, lname, &ctx.limit)); 671 PetscCheck(ctx.limit, PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "Limiter '%s' not found", lname); 672 673 /* Choose the physics from the list of registered models */ 674 { 675 PetscErrorCode (*r)(FVCtx *); 676 PetscCall(PetscFunctionListFind(physics, physname, &r)); 677 PetscCheck(r, PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "Physics '%s' not found", physname); 678 /* Create the physics, will set the number of fields and their names */ 679 PetscCall((*r)(&ctx)); 680 } 681 682 /* Create a DMDA to manage the parallel grid */ 683 PetscCall(DMDACreate1d(comm, DM_BOUNDARY_PERIODIC, 50, ctx.physics.dof, 2, NULL, &da)); 684 PetscCall(DMSetFromOptions(da)); 685 PetscCall(DMSetUp(da)); 686 /* Inform the DMDA of the field names provided by the physics. */ 687 /* The names will be shown in the title bars when run with -ts_monitor_draw_solution */ 688 for (i = 0; i < ctx.physics.dof; i++) PetscCall(DMDASetFieldName(da, i, ctx.physics.fieldname[i])); 689 PetscCall(DMDAGetInfo(da, 0, &Mx, 0, 0, 0, 0, 0, &dof, 0, 0, 0, 0, 0)); 690 PetscCall(DMDAGetCorners(da, &xs, 0, 0, &xm, 0, 0)); 691 692 /* Set coordinates of cell centers */ 693 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)); 694 695 /* Allocate work space for the Finite Volume solver (so it doesn't have to be reallocated on each function evaluation) */ 696 PetscCall(PetscMalloc4(dof * dof, &ctx.R, dof * dof, &ctx.Rinv, 2 * dof, &ctx.cjmpLR, 1 * dof, &ctx.cslope)); 697 PetscCall(PetscMalloc3(2 * dof, &ctx.uLR, dof, &ctx.flux, dof, &ctx.speeds)); 698 699 /* Create a vector to store the solution and to save the initial state */ 700 PetscCall(DMCreateGlobalVector(da, &X)); 701 PetscCall(VecDuplicate(X, &X0)); 702 PetscCall(VecDuplicate(X, &R)); 703 704 /*-------------------------------- create index for slow parts and fast parts ----------------------------------------*/ 705 PetscCall(PetscMalloc1(xm * dof, &index_slow)); 706 PetscCall(PetscMalloc1(xm * dof, &index_fast)); 707 for (i = xs; i < xs + xm; i++) { 708 if (ctx.xmin + i * (ctx.xmax - ctx.xmin) / (PetscReal)Mx + 0.5 * (ctx.xmax - ctx.xmin) / (PetscReal)Mx < 0) 709 for (k = 0; k < dof; k++) index_slow[islow++] = i * dof + k; 710 else 711 for (k = 0; k < dof; k++) index_fast[ifast++] = i * dof + k; 712 } /* this step need to be changed based on discontinuous point of a */ 713 PetscCall(ISCreateGeneral(PETSC_COMM_WORLD, islow, index_slow, PETSC_COPY_VALUES, &ctx.iss)); 714 PetscCall(ISCreateGeneral(PETSC_COMM_WORLD, ifast, index_fast, PETSC_COPY_VALUES, &ctx.isf)); 715 716 /* Create a time-stepping object */ 717 PetscCall(TSCreate(comm, &ts)); 718 PetscCall(TSSetDM(ts, da)); 719 PetscCall(TSSetRHSFunction(ts, R, FVRHSFunction, &ctx)); 720 PetscCall(TSRHSSplitSetIS(ts, "slow", ctx.iss)); 721 PetscCall(TSRHSSplitSetIS(ts, "fast", ctx.isf)); 722 PetscCall(TSRHSSplitSetRHSFunction(ts, "slow", NULL, FVRHSFunctionslow, &ctx)); 723 PetscCall(TSRHSSplitSetRHSFunction(ts, "fast", NULL, FVRHSFunctionfast, &ctx)); 724 725 if (ctx.recursive) { 726 TS subts; 727 islow = 0; 728 ifast = 0; 729 for (i = xs; i < xs + xm; i++) { 730 PetscReal coord = ctx.xmin + i * (ctx.xmax - ctx.xmin) / (PetscReal)Mx + 0.5 * (ctx.xmax - ctx.xmin) / (PetscReal)Mx; 731 if (coord >= 0 && coord < ctx.xmin + (ctx.xmax - ctx.xmin) * 3 / 4.) 732 for (k = 0; k < dof; k++) index_slow[islow++] = i * dof + k; 733 if (coord >= ctx.xmin + (ctx.xmax - ctx.xmin) * 3 / 4.) 734 for (k = 0; k < dof; k++) index_fast[ifast++] = i * dof + k; 735 } /* this step need to be changed based on discontinuous point of a */ 736 PetscCall(ISCreateGeneral(PETSC_COMM_WORLD, islow, index_slow, PETSC_COPY_VALUES, &ctx.iss2)); 737 PetscCall(ISCreateGeneral(PETSC_COMM_WORLD, ifast, index_fast, PETSC_COPY_VALUES, &ctx.isf2)); 738 739 PetscCall(TSRHSSplitGetSubTS(ts, "fast", &subts)); 740 PetscCall(TSRHSSplitSetIS(subts, "slow", ctx.iss2)); 741 PetscCall(TSRHSSplitSetIS(subts, "fast", ctx.isf2)); 742 PetscCall(TSRHSSplitSetRHSFunction(subts, "slow", NULL, FVRHSFunctionslow2, &ctx)); 743 PetscCall(TSRHSSplitSetRHSFunction(subts, "fast", NULL, FVRHSFunctionfast2, &ctx)); 744 } 745 746 /*PetscCall(TSSetType(ts,TSSSP));*/ 747 PetscCall(TSSetType(ts, TSMPRK)); 748 PetscCall(TSSetMaxTime(ts, 10)); 749 PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER)); 750 751 /* Compute initial conditions and starting time step */ 752 PetscCall(FVSample(&ctx, da, 0, X0)); 753 PetscCall(FVRHSFunction(ts, 0, X0, X, (void *)&ctx)); /* Initial function evaluation, only used to determine max speed */ 754 PetscCall(VecCopy(X0, X)); /* The function value was not used so we set X=X0 again */ 755 PetscCall(TSSetTimeStep(ts, ctx.cfl / ctx.cfl_idt)); 756 PetscCall(TSSetFromOptions(ts)); /* Take runtime options */ 757 PetscCall(SolutionStatsView(da, X, PETSC_VIEWER_STDOUT_WORLD)); 758 { 759 PetscInt steps; 760 PetscScalar mass_initial, mass_final, mass_difference; 761 762 PetscCall(TSSolve(ts, X)); 763 PetscCall(TSGetSolveTime(ts, &ptime)); 764 PetscCall(TSGetStepNumber(ts, &steps)); 765 PetscCall(PetscPrintf(comm, "Final time %g, steps %" PetscInt_FMT "\n", (double)ptime, steps)); 766 /* calculate the total mass at initial time and final time */ 767 mass_initial = 0.0; 768 mass_final = 0.0; 769 PetscCall(VecSum(X0, &mass_initial)); 770 PetscCall(VecSum(X, &mass_final)); 771 mass_difference = (ctx.xmax - ctx.xmin) / (PetscScalar)Mx * (mass_final - mass_initial); 772 PetscCall(PetscPrintf(comm, "Mass difference %g\n", (double)mass_difference)); 773 if (ctx.simulation) { 774 PetscViewer fd; 775 char filename[PETSC_MAX_PATH_LEN] = "binaryoutput"; 776 Vec XR; 777 PetscReal nrm1, nrmsup; 778 PetscBool flg; 779 780 PetscCall(PetscOptionsGetString(NULL, NULL, "-f", filename, sizeof(filename), &flg)); 781 PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_USER, "Must indicate binary file with the -f option"); 782 PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, filename, FILE_MODE_READ, &fd)); 783 PetscCall(VecDuplicate(X0, &XR)); 784 PetscCall(VecLoad(XR, fd)); 785 PetscCall(PetscViewerDestroy(&fd)); 786 PetscCall(VecAYPX(XR, -1, X)); 787 PetscCall(VecNorm(XR, NORM_1, &nrm1)); 788 PetscCall(VecNorm(XR, NORM_INFINITY, &nrmsup)); 789 nrm1 /= Mx; 790 PetscCall(PetscPrintf(comm, "Error ||x-x_e||_1 %g ||x-x_e||_sup %g\n", (double)nrm1, (double)nrmsup)); 791 PetscCall(VecDestroy(&XR)); 792 } 793 } 794 795 PetscCall(SolutionStatsView(da, X, PETSC_VIEWER_STDOUT_WORLD)); 796 if (draw & 0x1) PetscCall(VecView(X0, PETSC_VIEWER_DRAW_WORLD)); 797 if (draw & 0x2) PetscCall(VecView(X, PETSC_VIEWER_DRAW_WORLD)); 798 if (draw & 0x4) { 799 Vec Y; 800 PetscCall(VecDuplicate(X, &Y)); 801 PetscCall(FVSample(&ctx, da, ptime, Y)); 802 PetscCall(VecAYPX(Y, -1, X)); 803 PetscCall(VecView(Y, PETSC_VIEWER_DRAW_WORLD)); 804 PetscCall(VecDestroy(&Y)); 805 } 806 807 if (view_final) { 808 PetscViewer viewer; 809 PetscCall(PetscViewerASCIIOpen(PETSC_COMM_WORLD, final_fname, &viewer)); 810 PetscCall(PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_MATLAB)); 811 PetscCall(VecView(X, viewer)); 812 PetscCall(PetscViewerPopFormat(viewer)); 813 PetscCall(PetscViewerDestroy(&viewer)); 814 } 815 816 /* Clean up */ 817 PetscCall((*ctx.physics.destroy)(ctx.physics.user)); 818 for (i = 0; i < ctx.physics.dof; i++) PetscCall(PetscFree(ctx.physics.fieldname[i])); 819 PetscCall(PetscFree4(ctx.R, ctx.Rinv, ctx.cjmpLR, ctx.cslope)); 820 PetscCall(PetscFree3(ctx.uLR, ctx.flux, ctx.speeds)); 821 PetscCall(ISDestroy(&ctx.iss)); 822 PetscCall(ISDestroy(&ctx.isf)); 823 PetscCall(ISDestroy(&ctx.iss2)); 824 PetscCall(ISDestroy(&ctx.isf2)); 825 PetscCall(VecDestroy(&X)); 826 PetscCall(VecDestroy(&X0)); 827 PetscCall(VecDestroy(&R)); 828 PetscCall(DMDestroy(&da)); 829 PetscCall(TSDestroy(&ts)); 830 PetscCall(PetscFree(index_slow)); 831 PetscCall(PetscFree(index_fast)); 832 PetscCall(PetscFunctionListDestroy(&limiters)); 833 PetscCall(PetscFunctionListDestroy(&physics)); 834 PetscCall(PetscFinalize()); 835 return 0; 836 } 837 838 /*TEST 839 build: 840 requires: !complex 841 depends: finitevolume1d.c 842 843 test: 844 args: -da_grid_x 60 -initial 1 -xmin -1 -xmax 1 -limit mc -physics_advect_a1 1 -physics_advect_a2 2 -ts_time_step 0.025 -ts_max_steps 24 -ts_type rk -ts_rk_type 2a -ts_rk_dtratio 2 -ts_rk_multirate -ts_use_splitrhsfunction 0 845 846 test: 847 suffix: 2 848 args: -da_grid_x 60 -initial 1 -xmin -1 -xmax 1 -limit mc -physics_advect_a1 1 -physics_advect_a2 2 -ts_time_step 0.025 -ts_max_steps 24 -ts_type rk -ts_rk_type 2a -ts_rk_dtratio 2 -ts_rk_multirate -ts_use_splitrhsfunction 1 849 output_file: output/ex5_1.out 850 851 test: 852 suffix: 3 853 args: -da_grid_x 60 -initial 1 -xmin -1 -xmax 1 -limit mc -physics_advect_a1 1 -physics_advect_a2 2 -ts_time_step 0.025 -ts_max_steps 24 -ts_type mprk -ts_mprk_type 2a22 -ts_use_splitrhsfunction 0 854 855 test: 856 suffix: 4 857 args: -da_grid_x 60 -initial 1 -xmin -1 -xmax 1 -limit mc -physics_advect_a1 1 -physics_advect_a2 2 -ts_time_step 0.025 -ts_max_steps 24 -ts_type mprk -ts_mprk_type 2a22 -ts_use_splitrhsfunction 1 858 output_file: output/ex5_3.out 859 860 test: 861 suffix: 5 862 args: -da_grid_x 60 -initial 1 -xmin -1 -xmax 1 -limit mc -physics_advect_a1 1 -physics_advect_a2 2 -ts_time_step 0.025 -ts_max_steps 24 -ts_type rk -ts_rk_type 2a -ts_rk_dtratio 2 -ts_rk_multirate -ts_use_splitrhsfunction 0 -recursive_split 863 864 test: 865 suffix: 6 866 args: -da_grid_x 60 -initial 1 -xmin -1 -xmax 1 -limit mc -physics_advect_a1 1 -physics_advect_a2 2 -ts_time_step 0.025 -ts_max_steps 24 -ts_type rk -ts_rk_type 2a -ts_rk_dtratio 2 -ts_rk_multirate -ts_use_splitrhsfunction 1 -recursive_split 867 TEST*/ 868