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