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