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
RangeMod(PetscReal a,PetscReal xmin,PetscReal xmax)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 ----------------------------------- */
PhysicsDestroy_SimpleFree(void * vctx)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
PhysicsFlux_Advect(void * vctx,const PetscScalar * u,PetscScalar * flux,PetscReal * maxspeed)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
PhysicsSample_Advect(void * vctx,PetscInt initial,FVBCType bctype,PetscReal xmin,PetscReal xmax,PetscReal t,PetscReal x,PetscReal * u)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
PhysicsCreate_Advect(FVCtx * ctx)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
FVRHSFunction(TS ts,PetscReal time,Vec X,Vec F,void * vctx)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 PetscCallMPI(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)(1 / (2 * ctx->cfl_idt))));
313 }
314 PetscCall(PetscFree4(r, min, alpha, gamma));
315 PetscFunctionReturn(PETSC_SUCCESS);
316 }
317
FVRHSFunctionslow(TS ts,PetscReal time,Vec X,Vec F,void * vctx)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
FVRHSFunctionfast(TS ts,PetscReal time,Vec X,Vec F,void * vctx)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
FVSample(FVCtx * ctx,DM da,PetscReal time,Vec U)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
SolutionStatsView(DM da,Vec X,PetscViewer viewer)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 isascii;
606
607 PetscFunctionBeginUser;
608 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
609 PetscCheck(isascii, PETSC_COMM_SELF, PETSC_ERR_SUP, "Viewer type not supported");
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 PetscCallMPI(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 PetscFunctionReturn(PETSC_SUCCESS);
630 }
631
SolutionErrorNorms(FVCtx * ctx,DM da,PetscReal t,Vec X,PetscReal * nrm1)632 static PetscErrorCode SolutionErrorNorms(FVCtx *ctx, DM da, PetscReal t, Vec X, PetscReal *nrm1)
633 {
634 Vec Y;
635 PetscInt i, Mx, count_slow = 0, count_fast = 0;
636 const PetscScalar *ptr_X, *ptr_Y;
637
638 PetscFunctionBeginUser;
639 PetscCall(VecGetSize(X, &Mx));
640 PetscCall(VecDuplicate(X, &Y));
641 PetscCall(FVSample(ctx, da, t, Y));
642 const PetscReal hs = (ctx->xmax - ctx->xmin) / 2.0 * (ctx->hratio + 1.0) / Mx;
643 const PetscReal hf = (ctx->xmax - ctx->xmin) / 2.0 * (1.0 + 1.0 / ctx->hratio) / Mx;
644 count_slow = (PetscReal)Mx / (1.0 + ctx->hratio);
645 count_fast = Mx - count_slow;
646 PetscCall(VecGetArrayRead(X, &ptr_X));
647 PetscCall(VecGetArrayRead(Y, &ptr_Y));
648 for (i = 0; i < Mx; i++) {
649 if (i < count_slow / 2 || i > count_slow / 2 + count_fast - 1) *nrm1 += hs * PetscAbs(ptr_X[i] - ptr_Y[i]);
650 else *nrm1 += hf * PetscAbs(ptr_X[i] - ptr_Y[i]);
651 }
652 PetscCall(VecRestoreArrayRead(X, &ptr_X));
653 PetscCall(VecRestoreArrayRead(Y, &ptr_Y));
654 PetscCall(VecDestroy(&Y));
655 PetscFunctionReturn(PETSC_SUCCESS);
656 }
657
main(int argc,char * argv[])658 int main(int argc, char *argv[])
659 {
660 char physname[256] = "advect", final_fname[256] = "solution.m";
661 PetscFunctionList physics = 0;
662 MPI_Comm comm;
663 TS ts;
664 DM da;
665 Vec X, X0, R;
666 FVCtx ctx;
667 PetscInt i, k, dof, xs, xm, Mx, draw = 0, count_slow, count_fast, islow = 0, ifast = 0, *index_slow, *index_fast;
668 PetscBool view_final = PETSC_FALSE;
669 PetscReal ptime;
670
671 PetscFunctionBeginUser;
672 PetscCall(PetscInitialize(&argc, &argv, 0, help));
673 comm = PETSC_COMM_WORLD;
674 PetscCall(PetscMemzero(&ctx, sizeof(ctx)));
675
676 /* Register physical models to be available on the command line */
677 PetscCall(PetscFunctionListAdd(&physics, "advect", PhysicsCreate_Advect));
678
679 ctx.comm = comm;
680 ctx.cfl = 0.9;
681 ctx.bctype = FVBC_PERIODIC;
682 ctx.xmin = -1.0;
683 ctx.xmax = 1.0;
684 PetscOptionsBegin(comm, NULL, "Finite Volume solver options", "");
685 PetscCall(PetscOptionsReal("-xmin", "X min", "", ctx.xmin, &ctx.xmin, NULL));
686 PetscCall(PetscOptionsReal("-xmax", "X max", "", ctx.xmax, &ctx.xmax, NULL));
687 PetscCall(PetscOptionsInt("-draw", "Draw solution vector, bitwise OR of (1=initial,2=final,4=final error)", "", draw, &draw, NULL));
688 PetscCall(PetscOptionsString("-view_final", "Write final solution in ASCII MATLAB format to given file name", "", final_fname, final_fname, sizeof(final_fname), &view_final));
689 PetscCall(PetscOptionsInt("-initial", "Initial condition (depends on the physics)", "", ctx.initial, &ctx.initial, NULL));
690 PetscCall(PetscOptionsBool("-exact", "Compare errors with exact solution", "", ctx.exact, &ctx.exact, NULL));
691 PetscCall(PetscOptionsBool("-simulation", "Compare errors with reference solution", "", ctx.simulation, &ctx.simulation, NULL));
692 PetscCall(PetscOptionsReal("-cfl", "CFL number to time step at", "", ctx.cfl, &ctx.cfl, NULL));
693 PetscCall(PetscOptionsEnum("-bc_type", "Boundary condition", "", FVBCTypes, (PetscEnum)ctx.bctype, (PetscEnum *)&ctx.bctype, NULL));
694 PetscCall(PetscOptionsInt("-hratio", "Spacing ratio", "", ctx.hratio, &ctx.hratio, NULL));
695 PetscOptionsEnd();
696
697 /* Choose the physics from the list of registered models */
698 {
699 PetscErrorCode (*r)(FVCtx *);
700 PetscCall(PetscFunctionListFind(physics, physname, &r));
701 PetscCheck(r, PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "Physics '%s' not found", physname);
702 /* Create the physics, will set the number of fields and their names */
703 PetscCall((*r)(&ctx));
704 }
705
706 /* Create a DMDA to manage the parallel grid */
707 PetscCall(DMDACreate1d(comm, DM_BOUNDARY_PERIODIC, 50, ctx.physics.dof, 2, NULL, &da));
708 PetscCall(DMSetFromOptions(da));
709 PetscCall(DMSetUp(da));
710 /* Inform the DMDA of the field names provided by the physics. */
711 /* The names will be shown in the title bars when run with -ts_monitor_draw_solution */
712 for (i = 0; i < ctx.physics.dof; i++) PetscCall(DMDASetFieldName(da, i, ctx.physics.fieldname[i]));
713 PetscCall(DMDAGetInfo(da, 0, &Mx, 0, 0, 0, 0, 0, &dof, 0, 0, 0, 0, 0));
714 PetscCall(DMDAGetCorners(da, &xs, 0, 0, &xm, 0, 0));
715
716 /* Set coordinates of cell centers */
717 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));
718
719 /* Allocate work space for the Finite Volume solver (so it doesn't have to be reallocated on each function evaluation) */
720 PetscCall(PetscMalloc3(dof, &ctx.u, dof, &ctx.flux, dof, &ctx.speeds));
721
722 /* Create a vector to store the solution and to save the initial state */
723 PetscCall(DMCreateGlobalVector(da, &X));
724 PetscCall(VecDuplicate(X, &X0));
725 PetscCall(VecDuplicate(X, &R));
726
727 /* create index for slow parts and fast parts*/
728 count_slow = Mx / (1 + ctx.hratio);
729 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");
730 count_fast = Mx - count_slow;
731 ctx.sf = count_slow / 2;
732 ctx.fs = ctx.sf + count_fast;
733 PetscCall(PetscMalloc1(xm * dof, &index_slow));
734 PetscCall(PetscMalloc1(xm * dof, &index_fast));
735 for (i = xs; i < xs + xm; i++) {
736 if (i < count_slow / 2 || i > count_slow / 2 + count_fast - 1)
737 for (k = 0; k < dof; k++) index_slow[islow++] = i * dof + k;
738 else
739 for (k = 0; k < dof; k++) index_fast[ifast++] = i * dof + k;
740 }
741 PetscCall(ISCreateGeneral(PETSC_COMM_WORLD, islow, index_slow, PETSC_COPY_VALUES, &ctx.iss));
742 PetscCall(ISCreateGeneral(PETSC_COMM_WORLD, ifast, index_fast, PETSC_COPY_VALUES, &ctx.isf));
743
744 /* Create a time-stepping object */
745 PetscCall(TSCreate(comm, &ts));
746 PetscCall(TSSetDM(ts, da));
747 PetscCall(TSSetRHSFunction(ts, R, FVRHSFunction, &ctx));
748 PetscCall(TSRHSSplitSetIS(ts, "slow", ctx.iss));
749 PetscCall(TSRHSSplitSetIS(ts, "fast", ctx.isf));
750 PetscCall(TSRHSSplitSetRHSFunction(ts, "slow", NULL, FVRHSFunctionslow, &ctx));
751 PetscCall(TSRHSSplitSetRHSFunction(ts, "fast", NULL, FVRHSFunctionfast, &ctx));
752
753 PetscCall(TSSetType(ts, TSMPRK));
754 PetscCall(TSSetMaxTime(ts, 10));
755 PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER));
756
757 /* Compute initial conditions and starting time step */
758 PetscCall(FVSample(&ctx, da, 0, X0));
759 PetscCall(FVRHSFunction(ts, 0, X0, X, (void *)&ctx)); /* Initial function evaluation, only used to determine max speed */
760 PetscCall(VecCopy(X0, X)); /* The function value was not used so we set X=X0 again */
761 PetscCall(TSSetTimeStep(ts, ctx.cfl / ctx.cfl_idt));
762 PetscCall(TSSetFromOptions(ts)); /* Take runtime options */
763 PetscCall(SolutionStatsView(da, X, PETSC_VIEWER_STDOUT_WORLD));
764 {
765 PetscInt steps;
766 PetscScalar mass_initial, mass_final, mass_difference, mass_differenceg;
767 const PetscScalar *ptr_X, *ptr_X0;
768 const PetscReal hs = (ctx.xmax - ctx.xmin) / 2.0 / count_slow;
769 const PetscReal hf = (ctx.xmax - ctx.xmin) / 2.0 / count_fast;
770 PetscCall(TSSolve(ts, X));
771 PetscCall(TSGetSolveTime(ts, &ptime));
772 PetscCall(TSGetStepNumber(ts, &steps));
773 /* calculate the total mass at initial time and final time */
774 mass_initial = 0.0;
775 mass_final = 0.0;
776 PetscCall(DMDAVecGetArrayRead(da, X0, (void *)&ptr_X0));
777 PetscCall(DMDAVecGetArrayRead(da, X, (void *)&ptr_X));
778 for (i = xs; i < xs + xm; i++) {
779 if (i < ctx.sf || i > ctx.fs - 1) {
780 for (k = 0; k < dof; k++) {
781 mass_initial = mass_initial + hs * ptr_X0[i * dof + k];
782 mass_final = mass_final + hs * ptr_X[i * dof + k];
783 }
784 } else {
785 for (k = 0; k < dof; k++) {
786 mass_initial = mass_initial + hf * ptr_X0[i * dof + k];
787 mass_final = mass_final + hf * ptr_X[i * dof + k];
788 }
789 }
790 }
791 PetscCall(DMDAVecRestoreArrayRead(da, X0, (void *)&ptr_X0));
792 PetscCall(DMDAVecRestoreArrayRead(da, X, (void *)&ptr_X));
793 mass_difference = mass_final - mass_initial;
794 PetscCallMPI(MPIU_Allreduce(&mass_difference, &mass_differenceg, 1, MPIU_SCALAR, MPIU_SUM, comm));
795 PetscCall(PetscPrintf(comm, "Mass difference %g\n", (double)mass_differenceg));
796 PetscCall(PetscPrintf(comm, "Final time %g, steps %" PetscInt_FMT "\n", (double)ptime, steps));
797 if (ctx.exact) {
798 PetscReal nrm1 = 0;
799 PetscCall(SolutionErrorNorms(&ctx, da, ptime, X, &nrm1));
800 PetscCall(PetscPrintf(comm, "Error ||x-x_e||_1 %g\n", (double)nrm1));
801 }
802 if (ctx.simulation) {
803 PetscReal nrm1 = 0;
804 PetscViewer fd;
805 char filename[PETSC_MAX_PATH_LEN] = "binaryoutput";
806 Vec XR;
807 PetscBool flg;
808 const PetscScalar *ptr_XR;
809 PetscCall(PetscOptionsGetString(NULL, NULL, "-f", filename, sizeof(filename), &flg));
810 PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_USER, "Must indicate binary file with the -f option");
811 PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, filename, FILE_MODE_READ, &fd));
812 PetscCall(VecDuplicate(X0, &XR));
813 PetscCall(VecLoad(XR, fd));
814 PetscCall(PetscViewerDestroy(&fd));
815 PetscCall(VecGetArrayRead(X, &ptr_X));
816 PetscCall(VecGetArrayRead(XR, &ptr_XR));
817 for (i = 0; i < Mx; i++) {
818 if (i < count_slow / 2 || i > count_slow / 2 + count_fast - 1) nrm1 = nrm1 + hs * PetscAbs(ptr_X[i] - ptr_XR[i]);
819 else nrm1 = nrm1 + hf * PetscAbs(ptr_X[i] - ptr_XR[i]);
820 }
821 PetscCall(VecRestoreArrayRead(X, &ptr_X));
822 PetscCall(VecRestoreArrayRead(XR, &ptr_XR));
823 PetscCall(PetscPrintf(comm, "Error ||x-x_e||_1 %g\n", (double)nrm1));
824 PetscCall(VecDestroy(&XR));
825 }
826 }
827
828 PetscCall(SolutionStatsView(da, X, PETSC_VIEWER_STDOUT_WORLD));
829 if (draw & 0x1) PetscCall(VecView(X0, PETSC_VIEWER_DRAW_WORLD));
830 if (draw & 0x2) PetscCall(VecView(X, PETSC_VIEWER_DRAW_WORLD));
831 if (draw & 0x4) {
832 Vec Y;
833 PetscCall(VecDuplicate(X, &Y));
834 PetscCall(FVSample(&ctx, da, ptime, Y));
835 PetscCall(VecAYPX(Y, -1, X));
836 PetscCall(VecView(Y, PETSC_VIEWER_DRAW_WORLD));
837 PetscCall(VecDestroy(&Y));
838 }
839
840 if (view_final) {
841 PetscViewer viewer;
842 PetscCall(PetscViewerASCIIOpen(PETSC_COMM_WORLD, final_fname, &viewer));
843 PetscCall(PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_MATLAB));
844 PetscCall(VecView(X, viewer));
845 PetscCall(PetscViewerPopFormat(viewer));
846 PetscCall(PetscViewerDestroy(&viewer));
847 }
848
849 /* Clean up */
850 PetscCall((*ctx.physics.destroy)(ctx.physics.user));
851 for (i = 0; i < ctx.physics.dof; i++) PetscCall(PetscFree(ctx.physics.fieldname[i]));
852 PetscCall(PetscFree3(ctx.u, ctx.flux, ctx.speeds));
853 PetscCall(ISDestroy(&ctx.iss));
854 PetscCall(ISDestroy(&ctx.isf));
855 PetscCall(VecDestroy(&X));
856 PetscCall(VecDestroy(&X0));
857 PetscCall(VecDestroy(&R));
858 PetscCall(DMDestroy(&da));
859 PetscCall(TSDestroy(&ts));
860 PetscCall(PetscFree(index_slow));
861 PetscCall(PetscFree(index_fast));
862 PetscCall(PetscFunctionListDestroy(&physics));
863 PetscCall(PetscFinalize());
864 return 0;
865 }
866
867 /*TEST
868
869 build:
870 requires: !complex
871
872 test:
873 args: -da_grid_x 60 -initial 7 -xmin -1 -xmax 1 -hratio 2 -ts_time_step 0.025 -ts_max_steps 24 -ts_type rk -ts_rk_type 2a -ts_rk_dtratio 2 -ts_rk_multirate -ts_use_splitrhsfunction 0
874
875 test:
876 suffix: 2
877 args: -da_grid_x 60 -initial 7 -xmin -1 -xmax 1 -hratio 2 -ts_time_step 0.025 -ts_max_steps 24 -ts_type rk -ts_rk_type 2a -ts_rk_dtratio 2 -ts_rk_multirate -ts_use_splitrhsfunction 1
878 output_file: output/ex7_1.out
879
880 test:
881 suffix: 3
882 args: -da_grid_x 60 -initial 7 -xmin -1 -xmax 1 -hratio 2 -ts_time_step 0.025 -ts_max_steps 24 -ts_type mprk -ts_mprk_type 2a22 -ts_use_splitrhsfunction 0
883
884 test:
885 suffix: 4
886 args: -da_grid_x 60 -initial 7 -xmin -1 -xmax 1 -hratio 2 -ts_time_step 0.025 -ts_max_steps 24 -ts_type mprk -ts_mprk_type 2a22 -ts_use_splitrhsfunction 1
887 output_file: output/ex7_3.out
888
889 test:
890 suffix: 5
891 nsize: 2
892 args: -da_grid_x 60 -initial 7 -xmin -1 -xmax 1 -hratio 2 -ts_time_step 0.025 -ts_max_steps 24 -ts_type mprk -ts_mprk_type 2a22 -ts_use_splitrhsfunction 1
893 output_file: output/ex7_3.out
894 TEST*/
895