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