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