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