1 static const char help[] = "1D periodic Finite Volume solver by a particular slope limiter 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, say\n" 5 " hxs = (xmax - xmin)/2.0*(hratio+1.0)/Mx, \n" 6 " hxf = (xmax - xmin)/2.0*(1.0+1.0/hratio)/Mx, \n" 7 " with x belongs to (xmin,xmax), the number of total mesh points is Mx and the ratio between the meshsize of coarse\n\n" 8 " grids and fine grids is hratio.\n" 9 " exact - Exact Riemann solver which usually needs to perform a Newton iteration to connect\n" 10 " the states across shocks and rarefactions\n" 11 " simulation - use reference solution which is generated by smaller time step size to be true solution,\n" 12 " also the reference solution should be generated by user and stored in a binary file.\n" 13 " characteristic - Limit the characteristic variables, this is usually preferred (default)\n" 14 "Several initial conditions can be chosen with -initial N\n\n" 15 "The problem size should be set with -da_grid_x M\n\n" 16 "This script choose the slope limiter by biased second-order upwind procedure which is proposed by Van Leer in 1994\n" 17 " u(x_(k+1/2),t) = u(x_k,t) + phi(x_(k+1/2),t)*(u(x_k,t)-u(x_(k-1),t)) \n" 18 " limiter phi(x_(k+1/2),t) = max(0,min(r(k+1/2),min(2,gamma(k+1/2)*r(k+1/2)+alpha(k+1/2)))) \n" 19 " r(k+1/2) = (u(x_(k+1))-u(x_k))/(u(x_k)-u(x_(k-1))) \n" 20 " alpha(k+1/2) = (h_k*h_(k+1))/(h_(k-1)+h_k)/(h_(k-1)+h_k+h_(k+1)) \n" 21 " gamma(k+1/2) = h_k*(h_(k-1)+h_k)/(h_k+h_(k+1))/(h_(k-1)+h_k+h_(k+1)) \n"; 22 23 #include <petscts.h> 24 #include <petscdm.h> 25 #include <petscdmda.h> 26 #include <petscdraw.h> 27 #include <petscmath.h> 28 29 static inline PetscReal RangeMod(PetscReal a, PetscReal xmin, PetscReal xmax) { 30 PetscReal range = xmax - xmin; 31 return xmin + PetscFmodReal(range + PetscFmodReal(a, range), range); 32 } 33 34 /* --------------------------------- Finite Volume data structures ----------------------------------- */ 35 36 typedef enum { 37 FVBC_PERIODIC, 38 FVBC_OUTFLOW 39 } FVBCType; 40 static const char *FVBCTypes[] = {"PERIODIC", "OUTFLOW", "FVBCType", "FVBC_", 0}; 41 42 typedef struct { 43 PetscErrorCode (*sample)(void *, PetscInt, FVBCType, PetscReal, PetscReal, PetscReal, PetscReal, PetscReal *); 44 PetscErrorCode (*flux)(void *, const PetscScalar *, PetscScalar *, PetscReal *); 45 PetscErrorCode (*destroy)(void *); 46 void *user; 47 PetscInt dof; 48 char *fieldname[16]; 49 } PhysicsCtx; 50 51 typedef struct { 52 PhysicsCtx physics; 53 MPI_Comm comm; 54 char prefix[256]; 55 56 /* Local work arrays */ 57 PetscScalar *flux; /* Flux across interface */ 58 PetscReal *speeds; /* Speeds of each wave */ 59 PetscScalar *u; /* value at face */ 60 61 PetscReal cfl_idt; /* Max allowable value of 1/Delta t */ 62 PetscReal cfl; 63 PetscReal xmin, xmax; 64 PetscInt initial; 65 PetscBool exact; 66 PetscBool simulation; 67 FVBCType bctype; 68 PetscInt hratio; /* hratio = hslow/hfast */ 69 IS isf, iss; 70 PetscInt sf, fs; /* slow-fast and fast-slow interfaces */ 71 } FVCtx; 72 73 /* --------------------------------- Physics ----------------------------------- */ 74 static PetscErrorCode PhysicsDestroy_SimpleFree(void *vctx) { 75 PetscFunctionBeginUser; 76 PetscCall(PetscFree(vctx)); 77 PetscFunctionReturn(0); 78 } 79 80 /* --------------------------------- Advection ----------------------------------- */ 81 typedef struct { 82 PetscReal a; /* advective velocity */ 83 } AdvectCtx; 84 85 static PetscErrorCode PhysicsFlux_Advect(void *vctx, const PetscScalar *u, PetscScalar *flux, PetscReal *maxspeed) { 86 AdvectCtx *ctx = (AdvectCtx *)vctx; 87 PetscReal speed; 88 89 PetscFunctionBeginUser; 90 speed = ctx->a; 91 flux[0] = speed * u[0]; 92 *maxspeed = speed; 93 PetscFunctionReturn(0); 94 } 95 96 static PetscErrorCode PhysicsSample_Advect(void *vctx, PetscInt initial, FVBCType bctype, PetscReal xmin, PetscReal xmax, PetscReal t, PetscReal x, PetscReal *u) { 97 AdvectCtx *ctx = (AdvectCtx *)vctx; 98 PetscReal a = ctx->a, x0; 99 100 PetscFunctionBeginUser; 101 switch (bctype) { 102 case FVBC_OUTFLOW: x0 = x - a * t; break; 103 case FVBC_PERIODIC: x0 = RangeMod(x - a * t, xmin, xmax); break; 104 default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "unknown BCType"); 105 } 106 switch (initial) { 107 case 0: u[0] = (x0 < 0) ? 1 : -1; break; 108 case 1: u[0] = (x0 < 0) ? -1 : 1; break; 109 case 2: u[0] = (0 < x0 && x0 < 1) ? 1 : 0; break; 110 case 3: u[0] = PetscSinReal(2 * PETSC_PI * x0); break; 111 case 4: u[0] = PetscAbs(x0); break; 112 case 5: u[0] = (x0 < 0 || x0 > 0.5) ? 0 : PetscSqr(PetscSinReal(2 * PETSC_PI * x0)); break; 113 case 6: u[0] = (x0 < 0) ? 0 : ((x0 < 1) ? x0 : ((x0 < 2) ? 2 - x0 : 0)); break; 114 case 7: u[0] = PetscPowReal(PetscSinReal(PETSC_PI * x0), 10.0); break; 115 default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "unknown initial condition"); 116 } 117 PetscFunctionReturn(0); 118 } 119 120 static PetscErrorCode PhysicsCreate_Advect(FVCtx *ctx) { 121 AdvectCtx *user; 122 123 PetscFunctionBeginUser; 124 PetscCall(PetscNew(&user)); 125 ctx->physics.sample = PhysicsSample_Advect; 126 ctx->physics.flux = PhysicsFlux_Advect; 127 ctx->physics.destroy = PhysicsDestroy_SimpleFree; 128 ctx->physics.user = user; 129 ctx->physics.dof = 1; 130 PetscCall(PetscStrallocpy("u", &ctx->physics.fieldname[0])); 131 user->a = 1; 132 PetscOptionsBegin(ctx->comm, ctx->prefix, "Options for advection", ""); 133 { PetscCall(PetscOptionsReal("-physics_advect_a", "Speed", "", user->a, &user->a, NULL)); } 134 PetscOptionsEnd(); 135 PetscFunctionReturn(0); 136 } 137 138 /* --------------------------------- Finite Volume Solver ----------------------------------- */ 139 140 static PetscErrorCode FVRHSFunction(TS ts, PetscReal time, Vec X, Vec F, void *vctx) { 141 FVCtx *ctx = (FVCtx *)vctx; 142 PetscInt i, j, Mx, dof, xs, xm, sf = ctx->sf, fs = ctx->fs; 143 PetscReal hf, hs, cfl_idt = 0; 144 PetscScalar *x, *f, *r, *min, *alpha, *gamma; 145 Vec Xloc; 146 DM da; 147 148 PetscFunctionBeginUser; 149 PetscCall(TSGetDM(ts, &da)); 150 PetscCall(DMGetLocalVector(da, &Xloc)); /* Xloc contains ghost points */ 151 PetscCall(DMDAGetInfo(da, 0, &Mx, 0, 0, 0, 0, 0, &dof, 0, 0, 0, 0, 0)); /* Mx is the number of center points */ 152 hs = (ctx->xmax - ctx->xmin) / 2.0 * (ctx->hratio + 1.0) / Mx; 153 hf = (ctx->xmax - ctx->xmin) / 2.0 * (1.0 + 1.0 / ctx->hratio) / Mx; 154 PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, Xloc)); /* X is solution vector which does not contain ghost points */ 155 PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, Xloc)); 156 PetscCall(VecZeroEntries(F)); /* F is the right hand side function corresponds to center points */ 157 PetscCall(DMDAVecGetArray(da, Xloc, &x)); 158 PetscCall(DMDAVecGetArray(da, F, &f)); 159 PetscCall(DMDAGetCorners(da, &xs, 0, 0, &xm, 0, 0)); 160 PetscCall(PetscMalloc4(dof, &r, dof, &min, dof, &alpha, dof, &gamma)); 161 162 if (ctx->bctype == FVBC_OUTFLOW) { 163 for (i = xs - 2; i < 0; i++) { 164 for (j = 0; j < dof; j++) x[i * dof + j] = x[j]; 165 } 166 for (i = Mx; i < xs + xm + 2; i++) { 167 for (j = 0; j < dof; j++) x[i * dof + j] = x[(xs + xm - 1) * dof + j]; 168 } 169 } 170 171 for (i = xs; i < xs + xm + 1; i++) { 172 PetscReal maxspeed; 173 PetscScalar *u; 174 if (i < sf || i > fs + 1) { 175 u = &ctx->u[0]; 176 alpha[0] = 1.0 / 6.0; 177 gamma[0] = 1.0 / 3.0; 178 for (j = 0; j < dof; j++) { 179 r[j] = (x[i * dof + j] - x[(i - 1) * dof + j]) / (x[(i - 1) * dof + j] - x[(i - 2) * dof + j]); 180 min[j] = PetscMin(r[j], 2.0); 181 u[j] = x[(i - 1) * dof + j] + PetscMax(0, PetscMin(min[j], alpha[0] + gamma[0] * r[j])) * (x[(i - 1) * dof + j] - x[(i - 2) * dof + j]); 182 } 183 PetscCall((*ctx->physics.flux)(ctx->physics.user, u, ctx->flux, &maxspeed)); 184 cfl_idt = PetscMax(cfl_idt, PetscAbsScalar(maxspeed / hs)); 185 if (i > xs) { 186 for (j = 0; j < dof; j++) f[(i - 1) * dof + j] -= ctx->flux[j] / hs; 187 } 188 if (i < xs + xm) { 189 for (j = 0; j < dof; j++) f[i * dof + j] += ctx->flux[j] / hs; 190 } 191 } else if (i == sf) { 192 u = &ctx->u[0]; 193 alpha[0] = hs * hf / (hs + hs) / (hs + hs + hf); 194 gamma[0] = hs * (hs + hs) / (hs + hf) / (hs + hs + hf); 195 for (j = 0; j < dof; j++) { 196 r[j] = (x[i * dof + j] - x[(i - 1) * dof + j]) / (x[(i - 1) * dof + j] - x[(i - 2) * dof + j]); 197 min[j] = PetscMin(r[j], 2.0); 198 u[j] = x[(i - 1) * dof + j] + PetscMax(0, PetscMin(min[j], alpha[0] + gamma[0] * r[j])) * (x[(i - 1) * dof + j] - x[(i - 2) * dof + j]); 199 } 200 PetscCall((*ctx->physics.flux)(ctx->physics.user, u, ctx->flux, &maxspeed)); 201 if (i > xs) { 202 for (j = 0; j < dof; j++) f[(i - 1) * dof + j] -= ctx->flux[j] / hs; 203 } 204 if (i < xs + xm) { 205 for (j = 0; j < dof; j++) f[i * dof + j] += ctx->flux[j] / hf; 206 } 207 } else if (i == sf + 1) { 208 u = &ctx->u[0]; 209 alpha[0] = hf * hf / (hs + hf) / (hs + hf + hf); 210 gamma[0] = hf * (hs + hf) / (hf + hf) / (hs + hf + hf); 211 for (j = 0; j < dof; j++) { 212 r[j] = (x[i * dof + j] - x[(i - 1) * dof + j]) / (x[(i - 1) * dof + j] - x[(i - 2) * dof + j]); 213 min[j] = PetscMin(r[j], 2.0); 214 u[j] = x[(i - 1) * dof + j] + PetscMax(0, PetscMin(min[j], alpha[0] + gamma[0] * r[j])) * (x[(i - 1) * dof + j] - x[(i - 2) * dof + j]); 215 } 216 PetscCall((*ctx->physics.flux)(ctx->physics.user, u, ctx->flux, &maxspeed)); 217 if (i > xs) { 218 for (j = 0; j < dof; j++) f[(i - 1) * dof + j] -= ctx->flux[j] / hf; 219 } 220 if (i < xs + xm) { 221 for (j = 0; j < dof; j++) f[i * dof + j] += ctx->flux[j] / hf; 222 } 223 } else if (i > sf + 1 && i < fs) { 224 u = &ctx->u[0]; 225 alpha[0] = 1.0 / 6.0; 226 gamma[0] = 1.0 / 3.0; 227 for (j = 0; j < dof; j++) { 228 r[j] = (x[i * dof + j] - x[(i - 1) * dof + j]) / (x[(i - 1) * dof + j] - x[(i - 2) * dof + j]); 229 min[j] = PetscMin(r[j], 2.0); 230 u[j] = x[(i - 1) * dof + j] + PetscMax(0, PetscMin(min[j], alpha[0] + gamma[0] * r[j])) * (x[(i - 1) * dof + j] - x[(i - 2) * dof + j]); 231 } 232 PetscCall((*ctx->physics.flux)(ctx->physics.user, u, ctx->flux, &maxspeed)); 233 if (i > xs) { 234 for (j = 0; j < dof; j++) f[(i - 1) * dof + j] -= ctx->flux[j] / hf; 235 } 236 if (i < xs + xm) { 237 for (j = 0; j < dof; j++) f[i * dof + j] += ctx->flux[j] / hf; 238 } 239 } else if (i == fs) { 240 u = &ctx->u[0]; 241 alpha[0] = hf * hs / (hf + hf) / (hf + hf + hs); 242 gamma[0] = hf * (hf + hf) / (hf + hs) / (hf + hf + hs); 243 for (j = 0; j < dof; j++) { 244 r[j] = (x[i * dof + j] - x[(i - 1) * dof + j]) / (x[(i - 1) * dof + j] - x[(i - 2) * dof + j]); 245 min[j] = PetscMin(r[j], 2.0); 246 u[j] = x[(i - 1) * dof + j] + PetscMax(0, PetscMin(min[j], alpha[0] + gamma[0] * r[j])) * (x[(i - 1) * dof + j] - x[(i - 2) * dof + j]); 247 } 248 PetscCall((*ctx->physics.flux)(ctx->physics.user, u, ctx->flux, &maxspeed)); 249 if (i > xs) { 250 for (j = 0; j < dof; j++) f[(i - 1) * dof + j] -= ctx->flux[j] / hf; 251 } 252 if (i < xs + xm) { 253 for (j = 0; j < dof; j++) f[i * dof + j] += ctx->flux[j] / hs; 254 } 255 } else if (i == fs + 1) { 256 u = &ctx->u[0]; 257 alpha[0] = hs * hs / (hf + hs) / (hf + hs + hs); 258 gamma[0] = hs * (hf + hs) / (hs + hs) / (hf + hs + hs); 259 for (j = 0; j < dof; j++) { 260 r[j] = (x[i * dof + j] - x[(i - 1) * dof + j]) / (x[(i - 1) * dof + j] - x[(i - 2) * dof + j]); 261 min[j] = PetscMin(r[j], 2.0); 262 u[j] = x[(i - 1) * dof + j] + PetscMax(0, PetscMin(min[j], alpha[0] + gamma[0] * r[j])) * (x[(i - 1) * dof + j] - x[(i - 2) * dof + j]); 263 } 264 PetscCall((*ctx->physics.flux)(ctx->physics.user, u, ctx->flux, &maxspeed)); 265 if (i > xs) { 266 for (j = 0; j < dof; j++) f[(i - 1) * dof + j] -= ctx->flux[j] / hs; 267 } 268 if (i < xs + xm) { 269 for (j = 0; j < dof; j++) f[i * dof + j] += ctx->flux[j] / hs; 270 } 271 } 272 } 273 PetscCall(DMDAVecRestoreArray(da, Xloc, &x)); 274 PetscCall(DMDAVecRestoreArray(da, F, &f)); 275 PetscCall(DMRestoreLocalVector(da, &Xloc)); 276 PetscCallMPI(MPI_Allreduce(&cfl_idt, &ctx->cfl_idt, 1, MPIU_SCALAR, MPIU_MAX, PetscObjectComm((PetscObject)da))); 277 if (0) { 278 /* We need a way to inform the TS of a CFL constraint, this is a debugging fragment */ 279 PetscReal dt, tnow; 280 PetscCall(TSGetTimeStep(ts, &dt)); 281 PetscCall(TSGetTime(ts, &tnow)); 282 if (dt > 0.5 / ctx->cfl_idt) { PetscCall(PetscPrintf(ctx->comm, "Stability constraint exceeded at t=%g, dt %g > %g\n", (double)tnow, (double)dt, (double)(0.5 / ctx->cfl_idt))); } 283 } 284 PetscCall(PetscFree4(r, min, alpha, gamma)); 285 PetscFunctionReturn(0); 286 } 287 288 static PetscErrorCode FVRHSFunctionslow(TS ts, PetscReal time, Vec X, Vec F, void *vctx) { 289 FVCtx *ctx = (FVCtx *)vctx; 290 PetscInt i, j, Mx, dof, xs, xm, islow = 0, sf = ctx->sf, fs = ctx->fs; 291 PetscReal hf, hs; 292 PetscScalar *x, *f, *r, *min, *alpha, *gamma; 293 Vec Xloc; 294 DM da; 295 296 PetscFunctionBeginUser; 297 PetscCall(TSGetDM(ts, &da)); 298 PetscCall(DMGetLocalVector(da, &Xloc)); /* Xloc contains ghost points */ 299 PetscCall(DMDAGetInfo(da, 0, &Mx, 0, 0, 0, 0, 0, &dof, 0, 0, 0, 0, 0)); /* Mx is the number of center points */ 300 hs = (ctx->xmax - ctx->xmin) / 2.0 * (ctx->hratio + 1.0) / Mx; 301 hf = (ctx->xmax - ctx->xmin) / 2.0 * (1.0 + 1.0 / ctx->hratio) / Mx; 302 PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, Xloc)); /* X is solution vector which does not contain ghost points */ 303 PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, Xloc)); 304 PetscCall(VecZeroEntries(F)); /* F is the right hand side function corresponds to center points */ 305 PetscCall(DMDAVecGetArray(da, Xloc, &x)); 306 PetscCall(VecGetArray(F, &f)); 307 PetscCall(DMDAGetCorners(da, &xs, 0, 0, &xm, 0, 0)); 308 PetscCall(PetscMalloc4(dof, &r, dof, &min, dof, &alpha, dof, &gamma)); 309 310 if (ctx->bctype == FVBC_OUTFLOW) { 311 for (i = xs - 2; i < 0; i++) { 312 for (j = 0; j < dof; j++) x[i * dof + j] = x[j]; 313 } 314 for (i = Mx; i < xs + xm + 2; i++) { 315 for (j = 0; j < dof; j++) x[i * dof + j] = x[(xs + xm - 1) * dof + j]; 316 } 317 } 318 319 for (i = xs; i < xs + xm + 1; i++) { 320 PetscReal maxspeed; 321 PetscScalar *u; 322 if (i < sf) { 323 u = &ctx->u[0]; 324 alpha[0] = 1.0 / 6.0; 325 gamma[0] = 1.0 / 3.0; 326 for (j = 0; j < dof; j++) { 327 r[j] = (x[i * dof + j] - x[(i - 1) * dof + j]) / (x[(i - 1) * dof + j] - x[(i - 2) * dof + j]); 328 min[j] = PetscMin(r[j], 2.0); 329 u[j] = x[(i - 1) * dof + j] + PetscMax(0, PetscMin(min[j], alpha[0] + gamma[0] * r[j])) * (x[(i - 1) * dof + j] - x[(i - 2) * dof + j]); 330 } 331 PetscCall((*ctx->physics.flux)(ctx->physics.user, u, ctx->flux, &maxspeed)); 332 if (i > xs) { 333 for (j = 0; j < dof; j++) f[(islow - 1) * dof + j] -= ctx->flux[j] / hs; 334 } 335 if (i < xs + xm) { 336 for (j = 0; j < dof; j++) f[islow * dof + j] += ctx->flux[j] / hs; 337 islow++; 338 } 339 } else if (i == sf) { 340 u = &ctx->u[0]; 341 alpha[0] = hs * hf / (hs + hs) / (hs + hs + hf); 342 gamma[0] = hs * (hs + hs) / (hs + hf) / (hs + hs + hf); 343 for (j = 0; j < dof; j++) { 344 r[j] = (x[i * dof + j] - x[(i - 1) * dof + j]) / (x[(i - 1) * dof + j] - x[(i - 2) * dof + j]); 345 min[j] = PetscMin(r[j], 2.0); 346 u[j] = x[(i - 1) * dof + j] + PetscMax(0, PetscMin(min[j], alpha[0] + gamma[0] * r[j])) * (x[(i - 1) * dof + j] - x[(i - 2) * dof + j]); 347 } 348 PetscCall((*ctx->physics.flux)(ctx->physics.user, u, ctx->flux, &maxspeed)); 349 if (i > xs) { 350 for (j = 0; j < dof; j++) f[(islow - 1) * dof + j] -= ctx->flux[j] / hs; 351 } 352 } else if (i == fs) { 353 u = &ctx->u[0]; 354 alpha[0] = hf * hs / (hf + hf) / (hf + hf + hs); 355 gamma[0] = hf * (hf + hf) / (hf + hs) / (hf + hf + hs); 356 for (j = 0; j < dof; j++) { 357 r[j] = (x[i * dof + j] - x[(i - 1) * dof + j]) / (x[(i - 1) * dof + j] - x[(i - 2) * dof + j]); 358 min[j] = PetscMin(r[j], 2.0); 359 u[j] = x[(i - 1) * dof + j] + PetscMax(0, PetscMin(min[j], alpha[0] + gamma[0] * r[j])) * (x[(i - 1) * dof + j] - x[(i - 2) * dof + j]); 360 } 361 PetscCall((*ctx->physics.flux)(ctx->physics.user, u, ctx->flux, &maxspeed)); 362 if (i < xs + xm) { 363 for (j = 0; j < dof; j++) f[islow * dof + j] += ctx->flux[j] / hs; 364 islow++; 365 } 366 } else if (i == fs + 1) { 367 u = &ctx->u[0]; 368 alpha[0] = hs * hs / (hf + hs) / (hf + hs + hs); 369 gamma[0] = hs * (hf + hs) / (hs + hs) / (hf + hs + hs); 370 for (j = 0; j < dof; j++) { 371 r[j] = (x[i * dof + j] - x[(i - 1) * dof + j]) / (x[(i - 1) * dof + j] - x[(i - 2) * dof + j]); 372 min[j] = PetscMin(r[j], 2.0); 373 u[j] = x[(i - 1) * dof + j] + PetscMax(0, PetscMin(min[j], alpha[0] + gamma[0] * r[j])) * (x[(i - 1) * dof + j] - x[(i - 2) * dof + j]); 374 } 375 PetscCall((*ctx->physics.flux)(ctx->physics.user, u, ctx->flux, &maxspeed)); 376 if (i > xs) { 377 for (j = 0; j < dof; j++) f[(islow - 1) * dof + j] -= ctx->flux[j] / hs; 378 } 379 if (i < xs + xm) { 380 for (j = 0; j < dof; j++) f[islow * dof + j] += ctx->flux[j] / hs; 381 islow++; 382 } 383 } else if (i > fs + 1) { 384 u = &ctx->u[0]; 385 alpha[0] = 1.0 / 6.0; 386 gamma[0] = 1.0 / 3.0; 387 for (j = 0; j < dof; j++) { 388 r[j] = (x[i * dof + j] - x[(i - 1) * dof + j]) / (x[(i - 1) * dof + j] - x[(i - 2) * dof + j]); 389 min[j] = PetscMin(r[j], 2.0); 390 u[j] = x[(i - 1) * dof + j] + PetscMax(0, PetscMin(min[j], alpha[0] + gamma[0] * r[j])) * (x[(i - 1) * dof + j] - x[(i - 2) * dof + j]); 391 } 392 PetscCall((*ctx->physics.flux)(ctx->physics.user, u, ctx->flux, &maxspeed)); 393 if (i > xs) { 394 for (j = 0; j < dof; j++) f[(islow - 1) * dof + j] -= ctx->flux[j] / hs; 395 } 396 if (i < xs + xm) { 397 for (j = 0; j < dof; j++) f[islow * dof + j] += ctx->flux[j] / hs; 398 islow++; 399 } 400 } 401 } 402 PetscCall(DMDAVecRestoreArray(da, Xloc, &x)); 403 PetscCall(VecRestoreArray(F, &f)); 404 PetscCall(DMRestoreLocalVector(da, &Xloc)); 405 PetscCall(PetscFree4(r, min, alpha, gamma)); 406 PetscFunctionReturn(0); 407 } 408 409 static PetscErrorCode FVRHSFunctionfast(TS ts, PetscReal time, Vec X, Vec F, void *vctx) { 410 FVCtx *ctx = (FVCtx *)vctx; 411 PetscInt i, j, Mx, dof, xs, xm, ifast = 0, sf = ctx->sf, fs = ctx->fs; 412 PetscReal hf, hs; 413 PetscScalar *x, *f, *r, *min, *alpha, *gamma; 414 Vec Xloc; 415 DM da; 416 417 PetscFunctionBeginUser; 418 PetscCall(TSGetDM(ts, &da)); 419 PetscCall(DMGetLocalVector(da, &Xloc)); /* Xloc contains ghost points */ 420 PetscCall(DMDAGetInfo(da, 0, &Mx, 0, 0, 0, 0, 0, &dof, 0, 0, 0, 0, 0)); /* Mx is the number of center points */ 421 hs = (ctx->xmax - ctx->xmin) / 2.0 * (ctx->hratio + 1.0) / Mx; 422 hf = (ctx->xmax - ctx->xmin) / 2.0 * (1.0 + 1.0 / ctx->hratio) / Mx; 423 PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, Xloc)); /* X is solution vector which does not contain ghost points */ 424 PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, Xloc)); 425 PetscCall(VecZeroEntries(F)); /* F is the right hand side function corresponds to center points */ 426 PetscCall(DMDAVecGetArray(da, Xloc, &x)); 427 PetscCall(VecGetArray(F, &f)); 428 PetscCall(DMDAGetCorners(da, &xs, 0, 0, &xm, 0, 0)); 429 PetscCall(PetscMalloc4(dof, &r, dof, &min, dof, &alpha, dof, &gamma)); 430 431 if (ctx->bctype == FVBC_OUTFLOW) { 432 for (i = xs - 2; i < 0; i++) { 433 for (j = 0; j < dof; j++) x[i * dof + j] = x[j]; 434 } 435 for (i = Mx; i < xs + xm + 2; i++) { 436 for (j = 0; j < dof; j++) x[i * dof + j] = x[(xs + xm - 1) * dof + j]; 437 } 438 } 439 440 for (i = xs; i < xs + xm + 1; i++) { 441 PetscReal maxspeed; 442 PetscScalar *u; 443 if (i == sf) { 444 u = &ctx->u[0]; 445 alpha[0] = hs * hf / (hs + hs) / (hs + hs + hf); 446 gamma[0] = hs * (hs + hs) / (hs + hf) / (hs + hs + hf); 447 for (j = 0; j < dof; j++) { 448 r[j] = (x[i * dof + j] - x[(i - 1) * dof + j]) / (x[(i - 1) * dof + j] - x[(i - 2) * dof + j]); 449 min[j] = PetscMin(r[j], 2.0); 450 u[j] = x[(i - 1) * dof + j] + PetscMax(0, PetscMin(min[j], alpha[0] + gamma[0] * r[j])) * (x[(i - 1) * dof + j] - x[(i - 2) * dof + j]); 451 } 452 PetscCall((*ctx->physics.flux)(ctx->physics.user, u, ctx->flux, &maxspeed)); 453 if (i < xs + xm) { 454 for (j = 0; j < dof; j++) f[ifast * dof + j] += ctx->flux[j] / hf; 455 ifast++; 456 } 457 } else if (i == sf + 1) { 458 u = &ctx->u[0]; 459 alpha[0] = hf * hf / (hs + hf) / (hs + hf + hf); 460 gamma[0] = hf * (hs + hf) / (hf + hf) / (hs + hf + hf); 461 for (j = 0; j < dof; j++) { 462 r[j] = (x[i * dof + j] - x[(i - 1) * dof + j]) / (x[(i - 1) * dof + j] - x[(i - 2) * dof + j]); 463 min[j] = PetscMin(r[j], 2.0); 464 u[j] = x[(i - 1) * dof + j] + PetscMax(0, PetscMin(min[j], alpha[0] + gamma[0] * r[j])) * (x[(i - 1) * dof + j] - x[(i - 2) * dof + j]); 465 } 466 PetscCall((*ctx->physics.flux)(ctx->physics.user, u, ctx->flux, &maxspeed)); 467 if (i > xs) { 468 for (j = 0; j < dof; j++) f[(ifast - 1) * dof + j] -= ctx->flux[j] / hf; 469 } 470 if (i < xs + xm) { 471 for (j = 0; j < dof; j++) f[ifast * dof + j] += ctx->flux[j] / hf; 472 ifast++; 473 } 474 } else if (i > sf + 1 && i < fs) { 475 u = &ctx->u[0]; 476 alpha[0] = 1.0 / 6.0; 477 gamma[0] = 1.0 / 3.0; 478 for (j = 0; j < dof; j++) { 479 r[j] = (x[i * dof + j] - x[(i - 1) * dof + j]) / (x[(i - 1) * dof + j] - x[(i - 2) * dof + j]); 480 min[j] = PetscMin(r[j], 2.0); 481 u[j] = x[(i - 1) * dof + j] + PetscMax(0, PetscMin(min[j], alpha[0] + gamma[0] * r[j])) * (x[(i - 1) * dof + j] - x[(i - 2) * dof + j]); 482 } 483 PetscCall((*ctx->physics.flux)(ctx->physics.user, u, ctx->flux, &maxspeed)); 484 if (i > xs) { 485 for (j = 0; j < dof; j++) f[(ifast - 1) * dof + j] -= ctx->flux[j] / hf; 486 } 487 if (i < xs + xm) { 488 for (j = 0; j < dof; j++) f[ifast * dof + j] += ctx->flux[j] / hf; 489 ifast++; 490 } 491 } else if (i == fs) { 492 u = &ctx->u[0]; 493 alpha[0] = hf * hs / (hf + hf) / (hf + hf + hs); 494 gamma[0] = hf * (hf + hf) / (hf + hs) / (hf + hf + hs); 495 for (j = 0; j < dof; j++) { 496 r[j] = (x[i * dof + j] - x[(i - 1) * dof + j]) / (x[(i - 1) * dof + j] - x[(i - 2) * dof + j]); 497 min[j] = PetscMin(r[j], 2.0); 498 u[j] = x[(i - 1) * dof + j] + PetscMax(0, PetscMin(min[j], alpha[0] + gamma[0] * r[j])) * (x[(i - 1) * dof + j] - x[(i - 2) * dof + j]); 499 } 500 PetscCall((*ctx->physics.flux)(ctx->physics.user, u, ctx->flux, &maxspeed)); 501 if (i > xs) { 502 for (j = 0; j < dof; j++) f[(ifast - 1) * dof + j] -= ctx->flux[j] / hf; 503 } 504 } 505 } 506 PetscCall(DMDAVecRestoreArray(da, Xloc, &x)); 507 PetscCall(VecRestoreArray(F, &f)); 508 PetscCall(DMRestoreLocalVector(da, &Xloc)); 509 PetscCall(PetscFree4(r, min, alpha, gamma)); 510 PetscFunctionReturn(0); 511 } 512 513 /* --------------------------------- Finite Volume Solver for slow components ----------------------------------- */ 514 515 PetscErrorCode FVSample(FVCtx *ctx, DM da, PetscReal time, Vec U) { 516 PetscScalar *u, *uj, xj, xi; 517 PetscInt i, j, k, dof, xs, xm, Mx, count_slow, count_fast; 518 const PetscInt N = 200; 519 520 PetscFunctionBeginUser; 521 PetscCheck(ctx->physics.sample, PETSC_COMM_SELF, PETSC_ERR_SUP, "Physics has not provided a sampling function"); 522 PetscCall(DMDAGetInfo(da, 0, &Mx, 0, 0, 0, 0, 0, &dof, 0, 0, 0, 0, 0)); 523 PetscCall(DMDAGetCorners(da, &xs, 0, 0, &xm, 0, 0)); 524 PetscCall(DMDAVecGetArray(da, U, &u)); 525 PetscCall(PetscMalloc1(dof, &uj)); 526 const PetscReal hs = (ctx->xmax - ctx->xmin) / 2.0 * (ctx->hratio + 1.0) / Mx; 527 const PetscReal hf = (ctx->xmax - ctx->xmin) / 2.0 * (1.0 + 1.0 / ctx->hratio) / Mx; 528 count_slow = Mx / (1 + ctx->hratio); 529 count_fast = Mx - count_slow; 530 for (i = xs; i < xs + xm; i++) { 531 if (i * hs + 0.5 * hs < (ctx->xmax - ctx->xmin) * 0.25) { 532 xi = ctx->xmin + 0.5 * hs + i * hs; 533 /* Integrate over cell i using trapezoid rule with N points. */ 534 for (k = 0; k < dof; k++) u[i * dof + k] = 0; 535 for (j = 0; j < N + 1; j++) { 536 xj = xi + hs * (j - N / 2) / (PetscReal)N; 537 PetscCall((*ctx->physics.sample)(ctx->physics.user, ctx->initial, ctx->bctype, ctx->xmin, ctx->xmax, time, xj, uj)); 538 for (k = 0; k < dof; k++) u[i * dof + k] += ((j == 0 || j == N) ? 0.5 : 1.0) * uj[k] / N; 539 } 540 } else if ((ctx->xmax - ctx->xmin) * 0.25 + (i - count_slow / 2) * hf + 0.5 * hf < (ctx->xmax - ctx->xmin) * 0.75) { 541 xi = ctx->xmin + (ctx->xmax - ctx->xmin) * 0.25 + 0.5 * hf + (i - count_slow / 2) * hf; 542 /* Integrate over cell i using trapezoid rule with N points. */ 543 for (k = 0; k < dof; k++) u[i * dof + k] = 0; 544 for (j = 0; j < N + 1; j++) { 545 xj = xi + hf * (j - N / 2) / (PetscReal)N; 546 PetscCall((*ctx->physics.sample)(ctx->physics.user, ctx->initial, ctx->bctype, ctx->xmin, ctx->xmax, time, xj, uj)); 547 for (k = 0; k < dof; k++) u[i * dof + k] += ((j == 0 || j == N) ? 0.5 : 1.0) * uj[k] / N; 548 } 549 } else { 550 xi = ctx->xmin + (ctx->xmax - ctx->xmin) * 0.75 + 0.5 * hs + (i - count_slow / 2 - count_fast) * hs; 551 /* Integrate over cell i using trapezoid rule with N points. */ 552 for (k = 0; k < dof; k++) u[i * dof + k] = 0; 553 for (j = 0; j < N + 1; j++) { 554 xj = xi + hs * (j - N / 2) / (PetscReal)N; 555 PetscCall((*ctx->physics.sample)(ctx->physics.user, ctx->initial, ctx->bctype, ctx->xmin, ctx->xmax, time, xj, uj)); 556 for (k = 0; k < dof; k++) u[i * dof + k] += ((j == 0 || j == N) ? 0.5 : 1.0) * uj[k] / N; 557 } 558 } 559 } 560 PetscCall(DMDAVecRestoreArray(da, U, &u)); 561 PetscCall(PetscFree(uj)); 562 PetscFunctionReturn(0); 563 } 564 565 static PetscErrorCode SolutionStatsView(DM da, Vec X, PetscViewer viewer) { 566 PetscReal xmin, xmax; 567 PetscScalar sum, tvsum, tvgsum; 568 const PetscScalar *x; 569 PetscInt imin, imax, Mx, i, j, xs, xm, dof; 570 Vec Xloc; 571 PetscBool iascii; 572 573 PetscFunctionBeginUser; 574 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii)); 575 if (iascii) { 576 /* PETSc lacks a function to compute total variation norm (difficult in multiple dimensions), we do it here */ 577 PetscCall(DMGetLocalVector(da, &Xloc)); 578 PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, Xloc)); 579 PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, Xloc)); 580 PetscCall(DMDAVecGetArrayRead(da, Xloc, (void *)&x)); 581 PetscCall(DMDAGetCorners(da, &xs, 0, 0, &xm, 0, 0)); 582 PetscCall(DMDAGetInfo(da, 0, &Mx, 0, 0, 0, 0, 0, &dof, 0, 0, 0, 0, 0)); 583 tvsum = 0; 584 for (i = xs; i < xs + xm; i++) { 585 for (j = 0; j < dof; j++) tvsum += PetscAbsScalar(x[i * dof + j] - x[(i - 1) * dof + j]); 586 } 587 PetscCallMPI(MPI_Allreduce(&tvsum, &tvgsum, 1, MPIU_SCALAR, MPIU_SUM, PetscObjectComm((PetscObject)da))); 588 PetscCall(DMDAVecRestoreArrayRead(da, Xloc, (void *)&x)); 589 PetscCall(DMRestoreLocalVector(da, &Xloc)); 590 591 PetscCall(VecMin(X, &imin, &xmin)); 592 PetscCall(VecMax(X, &imax, &xmax)); 593 PetscCall(VecSum(X, &sum)); 594 PetscCall(PetscViewerASCIIPrintf(viewer, "Solution range [%g,%g] with minimum at %" PetscInt_FMT ", mean %g, ||x||_TV %g\n", (double)xmin, (double)xmax, imin, (double)(sum / Mx), (double)(tvgsum / Mx))); 595 } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Viewer type not supported"); 596 PetscFunctionReturn(0); 597 } 598 599 static PetscErrorCode SolutionErrorNorms(FVCtx *ctx, DM da, PetscReal t, Vec X, PetscReal *nrm1) { 600 Vec Y; 601 PetscInt i, Mx, count_slow = 0, count_fast = 0; 602 const PetscScalar *ptr_X, *ptr_Y; 603 604 PetscFunctionBeginUser; 605 PetscCall(VecGetSize(X, &Mx)); 606 PetscCall(VecDuplicate(X, &Y)); 607 PetscCall(FVSample(ctx, da, t, Y)); 608 const PetscReal hs = (ctx->xmax - ctx->xmin) / 2.0 * (ctx->hratio + 1.0) / Mx; 609 const PetscReal hf = (ctx->xmax - ctx->xmin) / 2.0 * (1.0 + 1.0 / ctx->hratio) / Mx; 610 count_slow = (PetscReal)Mx / (1.0 + ctx->hratio); 611 count_fast = Mx - count_slow; 612 PetscCall(VecGetArrayRead(X, &ptr_X)); 613 PetscCall(VecGetArrayRead(Y, &ptr_Y)); 614 for (i = 0; i < Mx; i++) { 615 if (i < count_slow / 2 || i > count_slow / 2 + count_fast - 1) *nrm1 += hs * PetscAbs(ptr_X[i] - ptr_Y[i]); 616 else *nrm1 += hf * PetscAbs(ptr_X[i] - ptr_Y[i]); 617 } 618 PetscCall(VecRestoreArrayRead(X, &ptr_X)); 619 PetscCall(VecRestoreArrayRead(Y, &ptr_Y)); 620 PetscCall(VecDestroy(&Y)); 621 PetscFunctionReturn(0); 622 } 623 624 int main(int argc, char *argv[]) { 625 char physname[256] = "advect", final_fname[256] = "solution.m"; 626 PetscFunctionList physics = 0; 627 MPI_Comm comm; 628 TS ts; 629 DM da; 630 Vec X, X0, R; 631 FVCtx ctx; 632 PetscInt i, k, dof, xs, xm, Mx, draw = 0, count_slow, count_fast, islow = 0, ifast = 0, *index_slow, *index_fast; 633 PetscBool view_final = PETSC_FALSE; 634 PetscReal ptime; 635 636 PetscFunctionBeginUser; 637 PetscCall(PetscInitialize(&argc, &argv, 0, help)); 638 comm = PETSC_COMM_WORLD; 639 PetscCall(PetscMemzero(&ctx, sizeof(ctx))); 640 641 /* Register physical models to be available on the command line */ 642 PetscCall(PetscFunctionListAdd(&physics, "advect", PhysicsCreate_Advect)); 643 644 ctx.comm = comm; 645 ctx.cfl = 0.9; 646 ctx.bctype = FVBC_PERIODIC; 647 ctx.xmin = -1.0; 648 ctx.xmax = 1.0; 649 PetscOptionsBegin(comm, NULL, "Finite Volume solver options", ""); 650 PetscCall(PetscOptionsReal("-xmin", "X min", "", ctx.xmin, &ctx.xmin, NULL)); 651 PetscCall(PetscOptionsReal("-xmax", "X max", "", ctx.xmax, &ctx.xmax, NULL)); 652 PetscCall(PetscOptionsInt("-draw", "Draw solution vector, bitwise OR of (1=initial,2=final,4=final error)", "", draw, &draw, NULL)); 653 PetscCall(PetscOptionsString("-view_final", "Write final solution in ASCII MATLAB format to given file name", "", final_fname, final_fname, sizeof(final_fname), &view_final)); 654 PetscCall(PetscOptionsInt("-initial", "Initial condition (depends on the physics)", "", ctx.initial, &ctx.initial, NULL)); 655 PetscCall(PetscOptionsBool("-exact", "Compare errors with exact solution", "", ctx.exact, &ctx.exact, NULL)); 656 PetscCall(PetscOptionsBool("-simulation", "Compare errors with reference solution", "", ctx.simulation, &ctx.simulation, NULL)); 657 PetscCall(PetscOptionsReal("-cfl", "CFL number to time step at", "", ctx.cfl, &ctx.cfl, NULL)); 658 PetscCall(PetscOptionsEnum("-bc_type", "Boundary condition", "", FVBCTypes, (PetscEnum)ctx.bctype, (PetscEnum *)&ctx.bctype, NULL)); 659 PetscCall(PetscOptionsInt("-hratio", "Spacing ratio", "", ctx.hratio, &ctx.hratio, NULL)); 660 PetscOptionsEnd(); 661 662 /* Choose the physics from the list of registered models */ 663 { 664 PetscErrorCode (*r)(FVCtx *); 665 PetscCall(PetscFunctionListFind(physics, physname, &r)); 666 PetscCheck(r, PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "Physics '%s' not found", physname); 667 /* Create the physics, will set the number of fields and their names */ 668 PetscCall((*r)(&ctx)); 669 } 670 671 /* Create a DMDA to manage the parallel grid */ 672 PetscCall(DMDACreate1d(comm, DM_BOUNDARY_PERIODIC, 50, ctx.physics.dof, 2, NULL, &da)); 673 PetscCall(DMSetFromOptions(da)); 674 PetscCall(DMSetUp(da)); 675 /* Inform the DMDA of the field names provided by the physics. */ 676 /* The names will be shown in the title bars when run with -ts_monitor_draw_solution */ 677 for (i = 0; i < ctx.physics.dof; i++) { PetscCall(DMDASetFieldName(da, i, ctx.physics.fieldname[i])); } 678 PetscCall(DMDAGetInfo(da, 0, &Mx, 0, 0, 0, 0, 0, &dof, 0, 0, 0, 0, 0)); 679 PetscCall(DMDAGetCorners(da, &xs, 0, 0, &xm, 0, 0)); 680 681 /* Set coordinates of cell centers */ 682 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)); 683 684 /* Allocate work space for the Finite Volume solver (so it doesn't have to be reallocated on each function evaluation) */ 685 PetscCall(PetscMalloc3(dof, &ctx.u, dof, &ctx.flux, dof, &ctx.speeds)); 686 687 /* Create a vector to store the solution and to save the initial state */ 688 PetscCall(DMCreateGlobalVector(da, &X)); 689 PetscCall(VecDuplicate(X, &X0)); 690 PetscCall(VecDuplicate(X, &R)); 691 692 /* create index for slow parts and fast parts*/ 693 count_slow = Mx / (1 + ctx.hratio); 694 PetscCheck(count_slow % 2 == 0, PETSC_COMM_WORLD, PETSC_ERR_USER, "Please adjust grid size Mx (-da_grid_x) and hratio (-hratio) so that Mx/(1+hartio) is even"); 695 count_fast = Mx - count_slow; 696 ctx.sf = count_slow / 2; 697 ctx.fs = ctx.sf + count_fast; 698 PetscCall(PetscMalloc1(xm * dof, &index_slow)); 699 PetscCall(PetscMalloc1(xm * dof, &index_fast)); 700 for (i = xs; i < xs + xm; i++) { 701 if (i < count_slow / 2 || i > count_slow / 2 + count_fast - 1) 702 for (k = 0; k < dof; k++) index_slow[islow++] = i * dof + k; 703 else 704 for (k = 0; k < dof; k++) index_fast[ifast++] = i * dof + k; 705 } 706 PetscCall(ISCreateGeneral(PETSC_COMM_WORLD, islow, index_slow, PETSC_COPY_VALUES, &ctx.iss)); 707 PetscCall(ISCreateGeneral(PETSC_COMM_WORLD, ifast, index_fast, PETSC_COPY_VALUES, &ctx.isf)); 708 709 /* Create a time-stepping object */ 710 PetscCall(TSCreate(comm, &ts)); 711 PetscCall(TSSetDM(ts, da)); 712 PetscCall(TSSetRHSFunction(ts, R, FVRHSFunction, &ctx)); 713 PetscCall(TSRHSSplitSetIS(ts, "slow", ctx.iss)); 714 PetscCall(TSRHSSplitSetIS(ts, "fast", ctx.isf)); 715 PetscCall(TSRHSSplitSetRHSFunction(ts, "slow", NULL, FVRHSFunctionslow, &ctx)); 716 PetscCall(TSRHSSplitSetRHSFunction(ts, "fast", NULL, FVRHSFunctionfast, &ctx)); 717 718 PetscCall(TSSetType(ts, TSMPRK)); 719 PetscCall(TSSetMaxTime(ts, 10)); 720 PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER)); 721 722 /* Compute initial conditions and starting time step */ 723 PetscCall(FVSample(&ctx, da, 0, X0)); 724 PetscCall(FVRHSFunction(ts, 0, X0, X, (void *)&ctx)); /* Initial function evaluation, only used to determine max speed */ 725 PetscCall(VecCopy(X0, X)); /* The function value was not used so we set X=X0 again */ 726 PetscCall(TSSetTimeStep(ts, ctx.cfl / ctx.cfl_idt)); 727 PetscCall(TSSetFromOptions(ts)); /* Take runtime options */ 728 PetscCall(SolutionStatsView(da, X, PETSC_VIEWER_STDOUT_WORLD)); 729 { 730 PetscInt steps; 731 PetscScalar mass_initial, mass_final, mass_difference, mass_differenceg; 732 const PetscScalar *ptr_X, *ptr_X0; 733 const PetscReal hs = (ctx.xmax - ctx.xmin) / 2.0 / count_slow; 734 const PetscReal hf = (ctx.xmax - ctx.xmin) / 2.0 / count_fast; 735 PetscCall(TSSolve(ts, X)); 736 PetscCall(TSGetSolveTime(ts, &ptime)); 737 PetscCall(TSGetStepNumber(ts, &steps)); 738 /* calculate the total mass at initial time and final time */ 739 mass_initial = 0.0; 740 mass_final = 0.0; 741 PetscCall(DMDAVecGetArrayRead(da, X0, (void *)&ptr_X0)); 742 PetscCall(DMDAVecGetArrayRead(da, X, (void *)&ptr_X)); 743 for (i = xs; i < xs + xm; i++) { 744 if (i < ctx.sf || i > ctx.fs - 1) { 745 for (k = 0; k < dof; k++) { 746 mass_initial = mass_initial + hs * ptr_X0[i * dof + k]; 747 mass_final = mass_final + hs * ptr_X[i * dof + k]; 748 } 749 } else { 750 for (k = 0; k < dof; k++) { 751 mass_initial = mass_initial + hf * ptr_X0[i * dof + k]; 752 mass_final = mass_final + hf * ptr_X[i * dof + k]; 753 } 754 } 755 } 756 PetscCall(DMDAVecRestoreArrayRead(da, X0, (void *)&ptr_X0)); 757 PetscCall(DMDAVecRestoreArrayRead(da, X, (void *)&ptr_X)); 758 mass_difference = mass_final - mass_initial; 759 PetscCallMPI(MPI_Allreduce(&mass_difference, &mass_differenceg, 1, MPIU_SCALAR, MPIU_SUM, comm)); 760 PetscCall(PetscPrintf(comm, "Mass difference %g\n", (double)mass_differenceg)); 761 PetscCall(PetscPrintf(comm, "Final time %g, steps %" PetscInt_FMT "\n", (double)ptime, steps)); 762 if (ctx.exact) { 763 PetscReal nrm1 = 0; 764 PetscCall(SolutionErrorNorms(&ctx, da, ptime, X, &nrm1)); 765 PetscCall(PetscPrintf(comm, "Error ||x-x_e||_1 %g\n", (double)nrm1)); 766 } 767 if (ctx.simulation) { 768 PetscReal nrm1 = 0; 769 PetscViewer fd; 770 char filename[PETSC_MAX_PATH_LEN] = "binaryoutput"; 771 Vec XR; 772 PetscBool flg; 773 const PetscScalar *ptr_XR; 774 PetscCall(PetscOptionsGetString(NULL, NULL, "-f", filename, sizeof(filename), &flg)); 775 PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_USER, "Must indicate binary file with the -f option"); 776 PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, filename, FILE_MODE_READ, &fd)); 777 PetscCall(VecDuplicate(X0, &XR)); 778 PetscCall(VecLoad(XR, fd)); 779 PetscCall(PetscViewerDestroy(&fd)); 780 PetscCall(VecGetArrayRead(X, &ptr_X)); 781 PetscCall(VecGetArrayRead(XR, &ptr_XR)); 782 for (i = 0; i < Mx; i++) { 783 if (i < count_slow / 2 || i > count_slow / 2 + count_fast - 1) nrm1 = nrm1 + hs * PetscAbs(ptr_X[i] - ptr_XR[i]); 784 else nrm1 = nrm1 + hf * PetscAbs(ptr_X[i] - ptr_XR[i]); 785 } 786 PetscCall(VecRestoreArrayRead(X, &ptr_X)); 787 PetscCall(VecRestoreArrayRead(XR, &ptr_XR)); 788 PetscCall(PetscPrintf(comm, "Error ||x-x_e||_1 %g\n", (double)nrm1)); 789 PetscCall(VecDestroy(&XR)); 790 } 791 } 792 793 PetscCall(SolutionStatsView(da, X, PETSC_VIEWER_STDOUT_WORLD)); 794 if (draw & 0x1) PetscCall(VecView(X0, PETSC_VIEWER_DRAW_WORLD)); 795 if (draw & 0x2) PetscCall(VecView(X, PETSC_VIEWER_DRAW_WORLD)); 796 if (draw & 0x4) { 797 Vec Y; 798 PetscCall(VecDuplicate(X, &Y)); 799 PetscCall(FVSample(&ctx, da, ptime, Y)); 800 PetscCall(VecAYPX(Y, -1, X)); 801 PetscCall(VecView(Y, PETSC_VIEWER_DRAW_WORLD)); 802 PetscCall(VecDestroy(&Y)); 803 } 804 805 if (view_final) { 806 PetscViewer viewer; 807 PetscCall(PetscViewerASCIIOpen(PETSC_COMM_WORLD, final_fname, &viewer)); 808 PetscCall(PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_MATLAB)); 809 PetscCall(VecView(X, viewer)); 810 PetscCall(PetscViewerPopFormat(viewer)); 811 PetscCall(PetscViewerDestroy(&viewer)); 812 } 813 814 /* Clean up */ 815 PetscCall((*ctx.physics.destroy)(ctx.physics.user)); 816 for (i = 0; i < ctx.physics.dof; i++) PetscCall(PetscFree(ctx.physics.fieldname[i])); 817 PetscCall(PetscFree3(ctx.u, ctx.flux, ctx.speeds)); 818 PetscCall(ISDestroy(&ctx.iss)); 819 PetscCall(ISDestroy(&ctx.isf)); 820 PetscCall(VecDestroy(&X)); 821 PetscCall(VecDestroy(&X0)); 822 PetscCall(VecDestroy(&R)); 823 PetscCall(DMDestroy(&da)); 824 PetscCall(TSDestroy(&ts)); 825 PetscCall(PetscFree(index_slow)); 826 PetscCall(PetscFree(index_fast)); 827 PetscCall(PetscFunctionListDestroy(&physics)); 828 PetscCall(PetscFinalize()); 829 return 0; 830 } 831 832 /*TEST 833 834 build: 835 requires: !complex 836 837 test: 838 args: -da_grid_x 60 -initial 7 -xmin -1 -xmax 1 -hratio 2 -ts_dt 0.025 -ts_max_steps 24 -ts_type rk -ts_rk_type 2a -ts_rk_dtratio 2 -ts_rk_multirate -ts_use_splitrhsfunction 0 839 840 test: 841 suffix: 2 842 args: -da_grid_x 60 -initial 7 -xmin -1 -xmax 1 -hratio 2 -ts_dt 0.025 -ts_max_steps 24 -ts_type rk -ts_rk_type 2a -ts_rk_dtratio 2 -ts_rk_multirate -ts_use_splitrhsfunction 1 843 output_file: output/ex7_1.out 844 845 test: 846 suffix: 3 847 args: -da_grid_x 60 -initial 7 -xmin -1 -xmax 1 -hratio 2 -ts_dt 0.025 -ts_max_steps 24 -ts_type mprk -ts_mprk_type 2a22 -ts_use_splitrhsfunction 0 848 849 test: 850 suffix: 4 851 args: -da_grid_x 60 -initial 7 -xmin -1 -xmax 1 -hratio 2 -ts_dt 0.025 -ts_max_steps 24 -ts_type mprk -ts_mprk_type 2a22 -ts_use_splitrhsfunction 1 852 output_file: output/ex7_3.out 853 854 test: 855 suffix: 5 856 nsize: 2 857 args: -da_grid_x 60 -initial 7 -xmin -1 -xmax 1 -hratio 2 -ts_dt 0.025 -ts_max_steps 24 -ts_type mprk -ts_mprk_type 2a22 -ts_use_splitrhsfunction 1 858 output_file: output/ex7_3.out 859 TEST*/ 860