1 /*
2   Note:
3   To performance a convergence study, a reference solution can be generated with a small stepsize and stored into a binary file, e.g. -ts_type ssp -ts_dt 1e-5 -ts_max_steps 30000 -ts_view_solution binary:reference.bin
4   Errors can be computed in the following runs with -simulation -f reference.bin
5 
6   Multirate RK options:
7   -ts_rk_dtratio is the ratio between larger time step size and small time step size
8   -ts_rk_multirate_type has three choices: none (for single step RK)
9                                            combined (for multirate method and user just need to provide the same RHS with the single step RK)
10                                            partitioned (for multiraet method and user need to provide two RHS, one is for fast components and the orther is for slow components
11 */
12 
13 static const char help[] = "1D periodic Finite Volume solver in slope-limiter form with semidiscrete time stepping.\n"
14   " advection   - Variable coefficient scalar advection\n"
15   "                u_t       + (a*u)_x               = 0\n"
16   " a is a piecewise function with two values say a0 and a1 (both a0 and a1 are constant).\n"
17   " in this toy problem we assume a=a0 when x<0 and a=a1 when x>0 with a0<a1 which has the same discontinuous point with initial condition.\n"
18   " we don't provide an exact solution, so you need to generate reference solution to do convergnece staudy,\n"
19   " more precisely, you need firstly generate a reference to a binary file say file.bin, then on commend line,\n"
20   " you should type -simulation -f file.bin.\n"
21   " you can choose the number of grids by -da_grid_x.\n"
22   " you can choose the value of a by -physics_advect_a1 and -physics_advect_a2.\n";
23 
24 #include <petscts.h>
25 #include <petscdm.h>
26 #include <petscdmda.h>
27 #include <petscdraw.h>
28 #include <petsc/private/tsimpl.h>
29 
30 #include <petsc/private/kernels/blockinvert.h> /* For the Kernel_*_gets_* stuff for BAIJ */
31 
32 #include "finitevolume1d.h"
33 
RangeMod(PetscReal a,PetscReal xmin,PetscReal xmax)34 PETSC_STATIC_INLINE PetscReal RangeMod(PetscReal a,PetscReal xmin,PetscReal xmax) { PetscReal range = xmax-xmin; return xmin +PetscFmodReal(range+PetscFmodReal(a,range),range); }
35 
36 /* --------------------------------- Advection ----------------------------------- */
37 
38 typedef struct {
39   PetscReal a[2];                  /* advective velocity */
40 } AdvectCtx;
41 
PhysicsRiemann_Advect(void * vctx,PetscInt m,const PetscScalar * uL,const PetscScalar * uR,PetscScalar * flux,PetscReal * maxspeed,PetscReal x,PetscReal xmin,PetscReal xmax)42 static PetscErrorCode PhysicsRiemann_Advect(void *vctx,PetscInt m,const PetscScalar *uL,const PetscScalar *uR,PetscScalar *flux,PetscReal *maxspeed,PetscReal x,PetscReal xmin, PetscReal xmax)
43 {
44   AdvectCtx *ctx = (AdvectCtx*)vctx;
45   PetscReal *speed;
46 
47   PetscFunctionBeginUser;
48   speed = ctx->a;
49   if (x==0 || x == xmin || x == xmax) flux[0] = PetscMax(0,(speed[0]+speed[1])/2.0)*uL[0] + PetscMin(0,(speed[0]+speed[1])/2.0)*uR[0]; /* if condition need to be changed base on different problem, '0' is the discontinuous point of a */
50   else if (x<0) flux[0] = PetscMax(0,speed[0])*uL[0] + PetscMin(0,speed[0])*uR[0];  /* else if condition need to be changed based on diferent problem, 'x = 0' is discontinuous point of a */
51   else flux[0] = PetscMax(0,speed[1])*uL[0] + PetscMin(0,speed[1])*uR[0];
52   *maxspeed = *speed;
53   PetscFunctionReturn(0);
54 }
55 
PhysicsCharacteristic_Advect(void * vctx,PetscInt m,const PetscScalar * u,PetscScalar * X,PetscScalar * Xi,PetscReal * speeds,PetscReal x)56 static PetscErrorCode PhysicsCharacteristic_Advect(void *vctx,PetscInt m,const PetscScalar *u,PetscScalar *X,PetscScalar *Xi,PetscReal *speeds,PetscReal x)
57 {
58   AdvectCtx *ctx = (AdvectCtx*)vctx;
59 
60   PetscFunctionBeginUser;
61   X[0]      = 1.;
62   Xi[0]     = 1.;
63   if (x<0) speeds[0] = ctx->a[0];  /* x is discontinuous point of a */
64   else    speeds[0] = ctx->a[1];
65   PetscFunctionReturn(0);
66 }
67 
PhysicsSample_Advect(void * vctx,PetscInt initial,FVBCType bctype,PetscReal xmin,PetscReal xmax,PetscReal t,PetscReal x,PetscReal * u)68 static PetscErrorCode PhysicsSample_Advect(void *vctx,PetscInt initial,FVBCType bctype,PetscReal xmin,PetscReal xmax,PetscReal t,PetscReal x,PetscReal *u)
69 {
70   AdvectCtx *ctx = (AdvectCtx*)vctx;
71   PetscReal *a    = ctx->a,x0;
72 
73   PetscFunctionBeginUser;
74   if (x<0){   /* x is cell center */
75     switch (bctype) {
76       case FVBC_OUTFLOW:  x0 = x-a[0]*t; break;
77       case FVBC_PERIODIC: x0 = RangeMod(x-a[0]*t,xmin,xmax); break;
78       default: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"unknown BCType");
79     }
80     switch (initial) {
81       case 0: u[0] = (x0 < 0) ? 1 : -1; break;
82       case 1: u[0] = (x0 < 0) ? -1 : 1; break;
83       case 2: u[0] = (0 < x0 && x0 < 1) ? 1 : 0; break;
84       case 3: u[0] = PetscSinReal(2*PETSC_PI*x0); break;
85       case 4: u[0] = PetscAbs(x0); break;
86       case 5: u[0] = (x0 < 0 || x0 > 0.5) ? 0 : PetscSqr(PetscSinReal(2*PETSC_PI*x0)); break;
87       case 6: u[0] = (x0 < 0) ? 0 : ((x0 < 1) ? x0 : ((x0 < 2) ? 2-x0 : 0)); break;
88       case 7: u[0] = PetscPowReal(PetscSinReal(PETSC_PI*x0),10.0);break;
89       default: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"unknown initial condition");
90     }
91   }
92   else{
93     switch (bctype) {
94       case FVBC_OUTFLOW:  x0 = x-a[1]*t; break;
95       case FVBC_PERIODIC: x0 = RangeMod(x-a[1]*t,xmin,xmax); break;
96       default: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"unknown BCType");
97     }
98     switch (initial) {
99       case 0: u[0] = (x0 < 0) ? 1 : -1; break;
100       case 1: u[0] = (x0 < 0) ? -1 : 1; break;
101       case 2: u[0] = (0 < x0 && x0 < 1) ? 1 : 0; break;
102       case 3: u[0] = PetscSinReal(2*PETSC_PI*x0); break;
103       case 4: u[0] = PetscAbs(x0); break;
104       case 5: u[0] = (x0 < 0 || x0 > 0.5) ? 0 : PetscSqr(PetscSinReal(2*PETSC_PI*x0)); break;
105       case 6: u[0] = (x0 < 0) ? 0 : ((x0 < 1) ? x0 : ((x0 < 2) ? 2-x0 : 0)); break;
106       case 7: u[0] = PetscPowReal(PetscSinReal(PETSC_PI*x0),10.0);break;
107       default: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"unknown initial condition");
108     }
109   }
110   PetscFunctionReturn(0);
111 }
112 
PhysicsCreate_Advect(FVCtx * ctx)113 static PetscErrorCode PhysicsCreate_Advect(FVCtx *ctx)
114 {
115   PetscErrorCode ierr;
116   AdvectCtx      *user;
117 
118   PetscFunctionBeginUser;
119   ierr = PetscNew(&user);CHKERRQ(ierr);
120   ctx->physics.sample         = PhysicsSample_Advect;
121   ctx->physics.riemann        = PhysicsRiemann_Advect;
122   ctx->physics.characteristic = PhysicsCharacteristic_Advect;
123   ctx->physics.destroy        = PhysicsDestroy_SimpleFree;
124   ctx->physics.user           = user;
125   ctx->physics.dof            = 1;
126   ierr = PetscStrallocpy("u",&ctx->physics.fieldname[0]);CHKERRQ(ierr);
127   user->a[0] = 1;
128   user->a[1] = 1;
129   ierr = PetscOptionsBegin(ctx->comm,ctx->prefix,"Options for advection","");CHKERRQ(ierr);
130   {
131     ierr = PetscOptionsReal("-physics_advect_a1","Speed1","",user->a[0],&user->a[0],NULL);CHKERRQ(ierr);
132     ierr = PetscOptionsReal("-physics_advect_a2","Speed2","",user->a[1],&user->a[1],NULL);CHKERRQ(ierr);
133   }
134   ierr = PetscOptionsEnd();CHKERRQ(ierr);
135   PetscFunctionReturn(0);
136 }
137 
138 /* --------------------------------- Finite Volume Solver for slow parts ----------------------------------- */
139 
FVRHSFunctionslow(TS ts,PetscReal time,Vec X,Vec F,void * vctx)140 PetscErrorCode FVRHSFunctionslow(TS ts,PetscReal time,Vec X,Vec F,void *vctx)
141 {
142   FVCtx          *ctx = (FVCtx*)vctx;
143   PetscErrorCode ierr;
144   PetscInt       i,j,k,Mx,dof,xs,xm,len_slow;
145   PetscReal      hx,cfl_idt = 0;
146   PetscScalar    *x,*f,*slope;
147   Vec            Xloc;
148   DM             da;
149 
150   PetscFunctionBeginUser;
151   ierr = TSGetDM(ts,&da);CHKERRQ(ierr);
152   ierr = DMGetLocalVector(da,&Xloc);CHKERRQ(ierr);
153   ierr = DMDAGetInfo(da,0, &Mx,0,0, 0,0,0, &dof,0,0,0,0,0);CHKERRQ(ierr);
154   hx   = (ctx->xmax-ctx->xmin)/Mx;
155   ierr = DMGlobalToLocalBegin(da,X,INSERT_VALUES,Xloc);CHKERRQ(ierr);
156   ierr = DMGlobalToLocalEnd(da,X,INSERT_VALUES,Xloc);CHKERRQ(ierr);
157   ierr = VecZeroEntries(F);CHKERRQ(ierr);
158   ierr = DMDAVecGetArray(da,Xloc,&x);CHKERRQ(ierr);
159   ierr = VecGetArray(F,&f);CHKERRQ(ierr);
160   ierr = DMDAGetArray(da,PETSC_TRUE,&slope);CHKERRQ(ierr);
161   ierr = DMDAGetCorners(da,&xs,0,0,&xm,0,0);CHKERRQ(ierr);
162   ierr = ISGetSize(ctx->iss,&len_slow);CHKERRQ(ierr);
163 
164   if (ctx->bctype == FVBC_OUTFLOW) {
165     for (i=xs-2; i<0; i++) {
166       for (j=0; j<dof; j++) x[i*dof+j] = x[j];
167     }
168     for (i=Mx; i<xs+xm+2; i++) {
169       for (j=0; j<dof; j++) x[i*dof+j] = x[(xs+xm-1)*dof+j];
170     }
171   }
172   for (i=xs-1; i<xs+xm+1; i++) {
173     struct _LimitInfo info;
174     PetscScalar       *cjmpL,*cjmpR;
175     if (i < len_slow+1) {
176       /* Determine the right eigenvectors R, where A = R \Lambda R^{-1} */
177       ierr = (*ctx->physics.characteristic)(ctx->physics.user,dof,&x[i*dof],ctx->R,ctx->Rinv,ctx->speeds,ctx->xmin+hx*i);CHKERRQ(ierr);
178       /* Evaluate jumps across interfaces (i-1, i) and (i, i+1), put in characteristic basis */
179       ierr  = PetscArrayzero(ctx->cjmpLR,2*dof);CHKERRQ(ierr);
180       cjmpL = &ctx->cjmpLR[0];
181       cjmpR = &ctx->cjmpLR[dof];
182       for (j=0; j<dof; j++) {
183         PetscScalar jmpL,jmpR;
184         jmpL = x[(i+0)*dof+j]-x[(i-1)*dof+j];
185         jmpR = x[(i+1)*dof+j]-x[(i+0)*dof+j];
186         for (k=0; k<dof; k++) {
187           cjmpL[k] += ctx->Rinv[k+j*dof]*jmpL;
188           cjmpR[k] += ctx->Rinv[k+j*dof]*jmpR;
189         }
190       }
191       /* Apply limiter to the left and right characteristic jumps */
192       info.m  = dof;
193       info.hx = hx;
194       (*ctx->limit)(&info,cjmpL,cjmpR,ctx->cslope);
195       for (j=0; j<dof; j++) ctx->cslope[j] /= hx; /* rescale to a slope */
196       for (j=0; j<dof; j++) {
197         PetscScalar tmp = 0;
198         for (k=0; k<dof; k++) tmp += ctx->R[j+k*dof]*ctx->cslope[k];
199         slope[i*dof+j] = tmp;
200       }
201     }
202   }
203 
204   for (i=xs; i<xs+xm+1; i++) {
205     PetscReal   maxspeed;
206     PetscScalar *uL,*uR;
207     if (i < len_slow) { /* slow parts can be changed based on a */
208       uL = &ctx->uLR[0];
209       uR = &ctx->uLR[dof];
210       for (j=0; j<dof; j++) {
211         uL[j] = x[(i-1)*dof+j]+slope[(i-1)*dof+j]*hx/2;
212         uR[j] = x[(i-0)*dof+j]-slope[(i-0)*dof+j]*hx/2;
213       }
214       ierr    = (*ctx->physics.riemann)(ctx->physics.user,dof,uL,uR,ctx->flux,&maxspeed,ctx->xmin+hx*i,ctx->xmin,ctx->xmax);CHKERRQ(ierr);
215       cfl_idt = PetscMax(cfl_idt,PetscAbsScalar(maxspeed/hx)); /* Max allowable value of 1/Delta t */
216       if (i > xs) {
217         for (j=0; j<dof; j++) f[(i-1)*dof+j] -= ctx->flux[j]/hx;
218       }
219       if (i < xs+xm) {
220         for (j=0; j<dof; j++) f[i*dof+j] += ctx->flux[j]/hx;
221       }
222     }
223     if (i == len_slow) { /* slow parts can be changed based on a */
224       uL = &ctx->uLR[0];
225       uR = &ctx->uLR[dof];
226       for (j=0; j<dof; j++) {
227         uL[j] = x[(i-1)*dof+j]+slope[(i-1)*dof+j]*hx/2;
228         uR[j] = x[(i-0)*dof+j]-slope[(i-0)*dof+j]*hx/2;
229       }
230       ierr    = (*ctx->physics.riemann)(ctx->physics.user,dof,uL,uR,ctx->flux,&maxspeed,ctx->xmin+hx*i,ctx->xmin,ctx->xmax);CHKERRQ(ierr);
231       cfl_idt = PetscMax(cfl_idt,PetscAbsScalar(maxspeed/hx)); /* Max allowable value of 1/Delta t */
232       if (i > xs) {
233         for (j=0; j<dof; j++) f[(i-1)*dof+j] -= ctx->flux[j]/hx;
234       }
235     }
236   }
237   ierr = DMDAVecRestoreArray(da,Xloc,&x);CHKERRQ(ierr);
238   ierr = VecRestoreArray(F,&f);CHKERRQ(ierr);
239   ierr = DMDARestoreArray(da,PETSC_TRUE,&slope);CHKERRQ(ierr);
240   ierr = DMRestoreLocalVector(da,&Xloc);CHKERRQ(ierr);
241   PetscFunctionReturn(0);
242 }
243 
244 /* --------------------------------- Finite Volume Solver for fast  parts ----------------------------------- */
245 
FVRHSFunctionfast(TS ts,PetscReal time,Vec X,Vec F,void * vctx)246 PetscErrorCode FVRHSFunctionfast(TS ts,PetscReal time,Vec X,Vec F,void *vctx)
247 {
248   FVCtx          *ctx = (FVCtx*)vctx;
249   PetscErrorCode ierr;
250   PetscInt       i,j,k,Mx,dof,xs,xm,len_slow;
251   PetscReal      hx,cfl_idt = 0;
252   PetscScalar    *x,*f,*slope;
253   Vec            Xloc;
254   DM             da;
255 
256   PetscFunctionBeginUser;
257   ierr = TSGetDM(ts,&da);CHKERRQ(ierr);
258   ierr = DMGetLocalVector(da,&Xloc);CHKERRQ(ierr);
259   ierr = DMDAGetInfo(da,0, &Mx,0,0, 0,0,0, &dof,0,0,0,0,0);CHKERRQ(ierr);
260   hx   = (ctx->xmax-ctx->xmin)/Mx;
261   ierr = DMGlobalToLocalBegin(da,X,INSERT_VALUES,Xloc);CHKERRQ(ierr);
262   ierr = DMGlobalToLocalEnd(da,X,INSERT_VALUES,Xloc);CHKERRQ(ierr);
263   ierr = VecZeroEntries(F);CHKERRQ(ierr);
264   ierr = DMDAVecGetArray(da,Xloc,&x);CHKERRQ(ierr);
265   ierr = VecGetArray(F,&f);CHKERRQ(ierr);
266   ierr = DMDAGetArray(da,PETSC_TRUE,&slope);CHKERRQ(ierr);
267   ierr = DMDAGetCorners(da,&xs,0,0,&xm,0,0);CHKERRQ(ierr);
268   ierr = ISGetSize(ctx->iss,&len_slow);CHKERRQ(ierr);
269 
270   if (ctx->bctype == FVBC_OUTFLOW) {
271     for (i=xs-2; i<0; i++) {
272       for (j=0; j<dof; j++) x[i*dof+j] = x[j];
273     }
274     for (i=Mx; i<xs+xm+2; i++) {
275       for (j=0; j<dof; j++) x[i*dof+j] = x[(xs+xm-1)*dof+j];
276     }
277   }
278   for (i=xs-1; i<xs+xm+1; i++) {
279     struct _LimitInfo info;
280     PetscScalar       *cjmpL,*cjmpR;
281     if (i > len_slow-2) {
282       ierr = (*ctx->physics.characteristic)(ctx->physics.user,dof,&x[i*dof],ctx->R,ctx->Rinv,ctx->speeds,ctx->xmin+hx*i);CHKERRQ(ierr);
283       ierr  = PetscArrayzero(ctx->cjmpLR,2*dof);CHKERRQ(ierr);
284       cjmpL = &ctx->cjmpLR[0];
285       cjmpR = &ctx->cjmpLR[dof];
286       for (j=0; j<dof; j++) {
287         PetscScalar jmpL,jmpR;
288         jmpL = x[(i+0)*dof+j]-x[(i-1)*dof+j];
289         jmpR = x[(i+1)*dof+j]-x[(i+0)*dof+j];
290         for (k=0; k<dof; k++) {
291           cjmpL[k] += ctx->Rinv[k+j*dof]*jmpL;
292           cjmpR[k] += ctx->Rinv[k+j*dof]*jmpR;
293         }
294       }
295       /* Apply limiter to the left and right characteristic jumps */
296       info.m  = dof;
297       info.hx = hx;
298       (*ctx->limit)(&info,cjmpL,cjmpR,ctx->cslope);
299       for (j=0; j<dof; j++) ctx->cslope[j] /= hx; /* rescale to a slope */
300       for (j=0; j<dof; j++) {
301         PetscScalar tmp = 0;
302         for (k=0; k<dof; k++) tmp += ctx->R[j+k*dof]*ctx->cslope[k];
303         slope[i*dof+j] = tmp;
304       }
305     }
306   }
307 
308   for (i=xs; i<xs+xm+1; i++) {
309     PetscReal   maxspeed;
310     PetscScalar *uL,*uR;
311     if (i > len_slow) {
312       uL = &ctx->uLR[0];
313       uR = &ctx->uLR[dof];
314       for (j=0; j<dof; j++) {
315         uL[j] = x[(i-1)*dof+j]+slope[(i-1)*dof+j]*hx/2;
316         uR[j] = x[(i-0)*dof+j]-slope[(i-0)*dof+j]*hx/2;
317       }
318       ierr    = (*ctx->physics.riemann)(ctx->physics.user,dof,uL,uR,ctx->flux,&maxspeed,ctx->xmin+hx*i,ctx->xmin,ctx->xmax);CHKERRQ(ierr);
319       cfl_idt = PetscMax(cfl_idt,PetscAbsScalar(maxspeed/hx)); /* Max allowable value of 1/Delta t */
320       if (i > xs) {
321         for (j=0; j<dof; j++) f[(i-len_slow-1)*dof+j] -= ctx->flux[j]/hx;
322       }
323       if (i < xs+xm) {
324         for (j=0; j<dof; j++) f[(i-len_slow)*dof+j] += ctx->flux[j]/hx;
325       }
326     }
327     if (i == len_slow) {
328       uL = &ctx->uLR[0];
329       uR = &ctx->uLR[dof];
330       for (j=0; j<dof; j++) {
331         uL[j] = x[(i-1)*dof+j]+slope[(i-1)*dof+j]*hx/2;
332         uR[j] = x[(i-0)*dof+j]-slope[(i-0)*dof+j]*hx/2;
333       }
334       ierr    = (*ctx->physics.riemann)(ctx->physics.user,dof,uL,uR,ctx->flux,&maxspeed,ctx->xmin+hx*i,ctx->xmin,ctx->xmax);CHKERRQ(ierr);
335       cfl_idt = PetscMax(cfl_idt,PetscAbsScalar(maxspeed/hx)); /* Max allowable value of 1/Delta t */
336       if (i < xs+xm) {
337         for (j=0; j<dof; j++) f[(i-len_slow)*dof+j] += ctx->flux[j]/hx;
338       }
339     }
340   }
341   ierr = DMDAVecRestoreArray(da,Xloc,&x);CHKERRQ(ierr);
342   ierr = VecRestoreArray(F,&f);CHKERRQ(ierr);
343   ierr = DMDARestoreArray(da,PETSC_TRUE,&slope);CHKERRQ(ierr);
344   ierr = DMRestoreLocalVector(da,&Xloc);CHKERRQ(ierr);
345   PetscFunctionReturn(0);
346 }
347 
FVRHSFunctionslow2(TS ts,PetscReal time,Vec X,Vec F,void * vctx)348 PetscErrorCode FVRHSFunctionslow2(TS ts,PetscReal time,Vec X,Vec F,void *vctx)
349 {
350   FVCtx          *ctx = (FVCtx*)vctx;
351   PetscErrorCode ierr;
352   PetscInt       i,j,k,Mx,dof,xs,xm,len_slow1,len_slow2;
353   PetscReal      hx,cfl_idt = 0;
354   PetscScalar    *x,*f,*slope;
355   Vec            Xloc;
356   DM             da;
357 
358   PetscFunctionBeginUser;
359   ierr = TSGetDM(ts,&da);CHKERRQ(ierr);
360   ierr = DMGetLocalVector(da,&Xloc);CHKERRQ(ierr);
361   ierr = DMDAGetInfo(da,0, &Mx,0,0, 0,0,0, &dof,0,0,0,0,0);CHKERRQ(ierr);
362   hx   = (ctx->xmax-ctx->xmin)/Mx;
363   ierr = DMGlobalToLocalBegin(da,X,INSERT_VALUES,Xloc);CHKERRQ(ierr);
364   ierr = DMGlobalToLocalEnd(da,X,INSERT_VALUES,Xloc);CHKERRQ(ierr);
365   ierr = VecZeroEntries(F);CHKERRQ(ierr);
366   ierr = DMDAVecGetArray(da,Xloc,&x);CHKERRQ(ierr);
367   ierr = VecGetArray(F,&f);CHKERRQ(ierr);
368   ierr = DMDAGetArray(da,PETSC_TRUE,&slope);CHKERRQ(ierr);
369   ierr = DMDAGetCorners(da,&xs,0,0,&xm,0,0);CHKERRQ(ierr);
370   ierr = ISGetSize(ctx->iss,&len_slow1);CHKERRQ(ierr);
371   ierr = ISGetSize(ctx->iss2,&len_slow2);CHKERRQ(ierr);
372   if (ctx->bctype == FVBC_OUTFLOW) {
373     for (i=xs-2; i<0; i++) {
374       for (j=0; j<dof; j++) x[i*dof+j] = x[j];
375     }
376     for (i=Mx; i<xs+xm+2; i++) {
377       for (j=0; j<dof; j++) x[i*dof+j] = x[(xs+xm-1)*dof+j];
378     }
379   }
380   for (i=xs-1; i<xs+xm+1; i++) {
381     struct _LimitInfo info;
382     PetscScalar       *cjmpL,*cjmpR;
383     if (i < len_slow1+len_slow2+1 && i > len_slow1-2) {
384       ierr = (*ctx->physics.characteristic)(ctx->physics.user,dof,&x[i*dof],ctx->R,ctx->Rinv,ctx->speeds,ctx->xmin+hx*i);CHKERRQ(ierr);
385       ierr  = PetscArrayzero(ctx->cjmpLR,2*dof);CHKERRQ(ierr);
386       cjmpL = &ctx->cjmpLR[0];
387       cjmpR = &ctx->cjmpLR[dof];
388       for (j=0; j<dof; j++) {
389         PetscScalar jmpL,jmpR;
390         jmpL = x[(i+0)*dof+j]-x[(i-1)*dof+j];
391         jmpR = x[(i+1)*dof+j]-x[(i+0)*dof+j];
392         for (k=0; k<dof; k++) {
393           cjmpL[k] += ctx->Rinv[k+j*dof]*jmpL;
394           cjmpR[k] += ctx->Rinv[k+j*dof]*jmpR;
395         }
396       }
397       /* Apply limiter to the left and right characteristic jumps */
398       info.m  = dof;
399       info.hx = hx;
400       (*ctx->limit)(&info,cjmpL,cjmpR,ctx->cslope);
401       for (j=0; j<dof; j++) ctx->cslope[j] /= hx; /* rescale to a slope */
402       for (j=0; j<dof; j++) {
403         PetscScalar tmp = 0;
404         for (k=0; k<dof; k++) tmp += ctx->R[j+k*dof]*ctx->cslope[k];
405         slope[i*dof+j] = tmp;
406       }
407     }
408   }
409 
410   for (i=xs; i<xs+xm+1; i++) {
411     PetscReal   maxspeed;
412     PetscScalar *uL,*uR;
413     if (i < len_slow1+len_slow2 && i > len_slow1) { /* slow parts can be changed based on a */
414       uL = &ctx->uLR[0];
415       uR = &ctx->uLR[dof];
416       for (j=0; j<dof; j++) {
417         uL[j] = x[(i-1)*dof+j]+slope[(i-1)*dof+j]*hx/2;
418         uR[j] = x[(i-0)*dof+j]-slope[(i-0)*dof+j]*hx/2;
419       }
420       ierr    = (*ctx->physics.riemann)(ctx->physics.user,dof,uL,uR,ctx->flux,&maxspeed,ctx->xmin+hx*i,ctx->xmin,ctx->xmax);CHKERRQ(ierr);
421       cfl_idt = PetscMax(cfl_idt,PetscAbsScalar(maxspeed/hx)); /* Max allowable value of 1/Delta t */
422       if (i > xs) {
423         for (j=0; j<dof; j++) f[(i-len_slow1-1)*dof+j] -= ctx->flux[j]/hx;
424       }
425       if (i < xs+xm) {
426         for (j=0; j<dof; j++) f[(i-len_slow1)*dof+j] += ctx->flux[j]/hx;
427       }
428     }
429     if (i == len_slow1+len_slow2) { /* slow parts can be changed based on a */
430       uL = &ctx->uLR[0];
431       uR = &ctx->uLR[dof];
432       for (j=0; j<dof; j++) {
433         uL[j] = x[(i-1)*dof+j]+slope[(i-1)*dof+j]*hx/2;
434         uR[j] = x[(i-0)*dof+j]-slope[(i-0)*dof+j]*hx/2;
435       }
436       ierr    = (*ctx->physics.riemann)(ctx->physics.user,dof,uL,uR,ctx->flux,&maxspeed,ctx->xmin+hx*i,ctx->xmin,ctx->xmax);CHKERRQ(ierr);
437       cfl_idt = PetscMax(cfl_idt,PetscAbsScalar(maxspeed/hx)); /* Max allowable value of 1/Delta t */
438       if (i > xs) {
439         for (j=0; j<dof; j++) f[(i-len_slow1-1)*dof+j] -= ctx->flux[j]/hx;
440       }
441     }
442     if (i == len_slow1) { /* slow parts can be changed based on a */
443       uL = &ctx->uLR[0];
444       uR = &ctx->uLR[dof];
445       for (j=0; j<dof; j++) {
446         uL[j] = x[(i-1)*dof+j]+slope[(i-1)*dof+j]*hx/2;
447         uR[j] = x[(i-0)*dof+j]-slope[(i-0)*dof+j]*hx/2;
448       }
449       ierr    = (*ctx->physics.riemann)(ctx->physics.user,dof,uL,uR,ctx->flux,&maxspeed,ctx->xmin+hx*i,ctx->xmin,ctx->xmax);CHKERRQ(ierr);
450       cfl_idt = PetscMax(cfl_idt,PetscAbsScalar(maxspeed/hx)); /* Max allowable value of 1/Delta t */
451       if (i < xs+xm) {
452         for (j=0; j<dof; j++) f[(i-len_slow1)*dof+j] += ctx->flux[j]/hx;
453       }
454     }
455   }
456 
457   ierr = DMDAVecRestoreArray(da,Xloc,&x);CHKERRQ(ierr);
458   ierr = VecRestoreArray(F,&f);CHKERRQ(ierr);
459   ierr = DMDARestoreArray(da,PETSC_TRUE,&slope);CHKERRQ(ierr);
460   ierr = DMRestoreLocalVector(da,&Xloc);CHKERRQ(ierr);
461   PetscFunctionReturn(0);
462 }
463 
FVRHSFunctionfast2(TS ts,PetscReal time,Vec X,Vec F,void * vctx)464 PetscErrorCode FVRHSFunctionfast2(TS ts,PetscReal time,Vec X,Vec F,void *vctx)
465 {
466   FVCtx          *ctx = (FVCtx*)vctx;
467   PetscErrorCode ierr;
468   PetscInt       i,j,k,Mx,dof,xs,xm,len_slow1,len_slow2;
469   PetscReal      hx,cfl_idt = 0;
470   PetscScalar    *x,*f,*slope;
471   Vec            Xloc;
472   DM             da;
473 
474   PetscFunctionBeginUser;
475   ierr = TSGetDM(ts,&da);CHKERRQ(ierr);
476   ierr = DMGetLocalVector(da,&Xloc);CHKERRQ(ierr);
477   ierr = DMDAGetInfo(da,0, &Mx,0,0, 0,0,0, &dof,0,0,0,0,0);CHKERRQ(ierr);
478   hx   = (ctx->xmax-ctx->xmin)/Mx;
479   ierr = DMGlobalToLocalBegin(da,X,INSERT_VALUES,Xloc);CHKERRQ(ierr);
480   ierr = DMGlobalToLocalEnd(da,X,INSERT_VALUES,Xloc);CHKERRQ(ierr);
481   ierr = VecZeroEntries(F);CHKERRQ(ierr);
482   ierr = DMDAVecGetArray(da,Xloc,&x);CHKERRQ(ierr);
483   ierr = VecGetArray(F,&f);CHKERRQ(ierr);
484   ierr = DMDAGetArray(da,PETSC_TRUE,&slope);CHKERRQ(ierr);
485   ierr = DMDAGetCorners(da,&xs,0,0,&xm,0,0);CHKERRQ(ierr);
486   ierr = ISGetSize(ctx->iss,&len_slow1);CHKERRQ(ierr);
487   ierr = ISGetSize(ctx->iss2,&len_slow2);CHKERRQ(ierr);
488 
489   if (ctx->bctype == FVBC_OUTFLOW) {
490     for (i=xs-2; i<0; i++) {
491       for (j=0; j<dof; j++) x[i*dof+j] = x[j];
492     }
493     for (i=Mx; i<xs+xm+2; i++) {
494       for (j=0; j<dof; j++) x[i*dof+j] = x[(xs+xm-1)*dof+j];
495     }
496   }
497   for (i=xs-1; i<xs+xm+1; i++) {
498     struct _LimitInfo info;
499     PetscScalar       *cjmpL,*cjmpR;
500     if (i > len_slow1+len_slow2-2) {
501       ierr = (*ctx->physics.characteristic)(ctx->physics.user,dof,&x[i*dof],ctx->R,ctx->Rinv,ctx->speeds,ctx->xmin+hx*i);CHKERRQ(ierr);
502       ierr  = PetscArrayzero(ctx->cjmpLR,2*dof);CHKERRQ(ierr);
503       cjmpL = &ctx->cjmpLR[0];
504       cjmpR = &ctx->cjmpLR[dof];
505       for (j=0; j<dof; j++) {
506         PetscScalar jmpL,jmpR;
507         jmpL = x[(i+0)*dof+j]-x[(i-1)*dof+j];
508         jmpR = x[(i+1)*dof+j]-x[(i+0)*dof+j];
509         for (k=0; k<dof; k++) {
510           cjmpL[k] += ctx->Rinv[k+j*dof]*jmpL;
511           cjmpR[k] += ctx->Rinv[k+j*dof]*jmpR;
512         }
513       }
514       /* Apply limiter to the left and right characteristic jumps */
515       info.m  = dof;
516       info.hx = hx;
517       (*ctx->limit)(&info,cjmpL,cjmpR,ctx->cslope);
518       for (j=0; j<dof; j++) ctx->cslope[j] /= hx; /* rescale to a slope */
519       for (j=0; j<dof; j++) {
520         PetscScalar tmp = 0;
521         for (k=0; k<dof; k++) tmp += ctx->R[j+k*dof]*ctx->cslope[k];
522         slope[i*dof+j] = tmp;
523       }
524     }
525   }
526 
527   for (i=xs; i<xs+xm+1; i++) {
528     PetscReal   maxspeed;
529     PetscScalar *uL,*uR;
530     if (i > len_slow1+len_slow2) {
531       uL = &ctx->uLR[0];
532       uR = &ctx->uLR[dof];
533       for (j=0; j<dof; j++) {
534         uL[j] = x[(i-1)*dof+j]+slope[(i-1)*dof+j]*hx/2;
535         uR[j] = x[(i-0)*dof+j]-slope[(i-0)*dof+j]*hx/2;
536       }
537       ierr    = (*ctx->physics.riemann)(ctx->physics.user,dof,uL,uR,ctx->flux,&maxspeed,ctx->xmin+hx*i,ctx->xmin,ctx->xmax);CHKERRQ(ierr);
538       cfl_idt = PetscMax(cfl_idt,PetscAbsScalar(maxspeed/hx)); /* Max allowable value of 1/Delta t */
539       if (i > xs) {
540         for (j=0; j<dof; j++) f[(i-len_slow1-len_slow2-1)*dof+j] -= ctx->flux[j]/hx;
541       }
542       if (i < xs+xm) {
543         for (j=0; j<dof; j++) f[(i-len_slow1-len_slow2)*dof+j] += ctx->flux[j]/hx;
544       }
545     }
546     if (i == len_slow1+len_slow2) {
547       uL = &ctx->uLR[0];
548       uR = &ctx->uLR[dof];
549       for (j=0; j<dof; j++) {
550         uL[j] = x[(i-1)*dof+j]+slope[(i-1)*dof+j]*hx/2;
551         uR[j] = x[(i-0)*dof+j]-slope[(i-0)*dof+j]*hx/2;
552       }
553       ierr    = (*ctx->physics.riemann)(ctx->physics.user,dof,uL,uR,ctx->flux,&maxspeed,ctx->xmin+hx*i,ctx->xmin,ctx->xmax);CHKERRQ(ierr);
554       cfl_idt = PetscMax(cfl_idt,PetscAbsScalar(maxspeed/hx)); /* Max allowable value of 1/Delta t */
555       if (i < xs+xm) {
556         for (j=0; j<dof; j++) f[(i-len_slow1-len_slow2)*dof+j] += ctx->flux[j]/hx;
557       }
558     }
559   }
560   ierr = DMDAVecRestoreArray(da,Xloc,&x);CHKERRQ(ierr);
561   ierr = VecRestoreArray(F,&f);CHKERRQ(ierr);
562   ierr = DMDARestoreArray(da,PETSC_TRUE,&slope);CHKERRQ(ierr);
563   ierr = DMRestoreLocalVector(da,&Xloc);CHKERRQ(ierr);
564   PetscFunctionReturn(0);
565 }
566 
main(int argc,char * argv[])567 int main(int argc,char *argv[])
568 {
569   char              lname[256] = "mc",physname[256] = "advect",final_fname[256] = "solution.m";
570   PetscFunctionList limiters   = 0,physics = 0;
571   MPI_Comm          comm;
572   TS                ts;
573   DM                da;
574   Vec               X,X0,R;
575   FVCtx             ctx;
576   PetscInt          i,k,dof,xs,xm,Mx,draw = 0,*index_slow,*index_fast,islow = 0,ifast = 0;
577   PetscBool         view_final = PETSC_FALSE;
578   PetscReal         ptime;
579   PetscErrorCode    ierr;
580 
581   ierr = PetscInitialize(&argc,&argv,0,help);if (ierr) return ierr;
582   comm = PETSC_COMM_WORLD;
583   ierr = PetscMemzero(&ctx,sizeof(ctx));CHKERRQ(ierr);
584 
585   /* Register limiters to be available on the command line */
586   ierr = PetscFunctionListAdd(&limiters,"upwind"              ,Limit_Upwind);CHKERRQ(ierr);
587   ierr = PetscFunctionListAdd(&limiters,"lax-wendroff"        ,Limit_LaxWendroff);CHKERRQ(ierr);
588   ierr = PetscFunctionListAdd(&limiters,"beam-warming"        ,Limit_BeamWarming);CHKERRQ(ierr);
589   ierr = PetscFunctionListAdd(&limiters,"fromm"               ,Limit_Fromm);CHKERRQ(ierr);
590   ierr = PetscFunctionListAdd(&limiters,"minmod"              ,Limit_Minmod);CHKERRQ(ierr);
591   ierr = PetscFunctionListAdd(&limiters,"superbee"            ,Limit_Superbee);CHKERRQ(ierr);
592   ierr = PetscFunctionListAdd(&limiters,"mc"                  ,Limit_MC);CHKERRQ(ierr);
593   ierr = PetscFunctionListAdd(&limiters,"vanleer"             ,Limit_VanLeer);CHKERRQ(ierr);
594   ierr = PetscFunctionListAdd(&limiters,"vanalbada"           ,Limit_VanAlbada);CHKERRQ(ierr);
595   ierr = PetscFunctionListAdd(&limiters,"vanalbadatvd"        ,Limit_VanAlbadaTVD);CHKERRQ(ierr);
596   ierr = PetscFunctionListAdd(&limiters,"koren"               ,Limit_Koren);CHKERRQ(ierr);
597   ierr = PetscFunctionListAdd(&limiters,"korensym"            ,Limit_KorenSym);CHKERRQ(ierr);
598   ierr = PetscFunctionListAdd(&limiters,"koren3"              ,Limit_Koren3);CHKERRQ(ierr);
599   ierr = PetscFunctionListAdd(&limiters,"cada-torrilhon2"     ,Limit_CadaTorrilhon2);CHKERRQ(ierr);
600   ierr = PetscFunctionListAdd(&limiters,"cada-torrilhon3-r0p1",Limit_CadaTorrilhon3R0p1);CHKERRQ(ierr);
601   ierr = PetscFunctionListAdd(&limiters,"cada-torrilhon3-r1"  ,Limit_CadaTorrilhon3R1);CHKERRQ(ierr);
602   ierr = PetscFunctionListAdd(&limiters,"cada-torrilhon3-r10" ,Limit_CadaTorrilhon3R10);CHKERRQ(ierr);
603   ierr = PetscFunctionListAdd(&limiters,"cada-torrilhon3-r100",Limit_CadaTorrilhon3R100);CHKERRQ(ierr);
604 
605   /* Register physical models to be available on the command line */
606   ierr = PetscFunctionListAdd(&physics,"advect"          ,PhysicsCreate_Advect);CHKERRQ(ierr);
607 
608   ctx.comm = comm;
609   ctx.cfl  = 0.9;
610   ctx.bctype = FVBC_PERIODIC;
611   ctx.xmin = -1.0;
612   ctx.xmax = 1.0;
613   ierr = PetscOptionsBegin(comm,NULL,"Finite Volume solver options","");CHKERRQ(ierr);
614   ierr = PetscOptionsReal("-xmin","X min","",ctx.xmin,&ctx.xmin,NULL);CHKERRQ(ierr);
615   ierr = PetscOptionsReal("-xmax","X max","",ctx.xmax,&ctx.xmax,NULL);CHKERRQ(ierr);
616   ierr = PetscOptionsFList("-limit","Name of flux limiter to use","",limiters,lname,lname,sizeof(lname),NULL);CHKERRQ(ierr);
617   ierr = PetscOptionsInt("-draw","Draw solution vector, bitwise OR of (1=initial,2=final,4=final error)","",draw,&draw,NULL);CHKERRQ(ierr);
618   ierr = PetscOptionsString("-view_final","Write final solution in ASCII MATLAB format to given file name","",final_fname,final_fname,sizeof(final_fname),&view_final);CHKERRQ(ierr);
619   ierr = PetscOptionsInt("-initial","Initial condition (depends on the physics)","",ctx.initial,&ctx.initial,NULL);CHKERRQ(ierr);
620   ierr = PetscOptionsBool("-simulation","Compare errors with reference solution","",ctx.simulation,&ctx.simulation,NULL);CHKERRQ(ierr);
621   ierr = PetscOptionsReal("-cfl","CFL number to time step at","",ctx.cfl,&ctx.cfl,NULL);CHKERRQ(ierr);
622   ierr = PetscOptionsEnum("-bc_type","Boundary condition","",FVBCTypes,(PetscEnum)ctx.bctype,(PetscEnum*)&ctx.bctype,NULL);CHKERRQ(ierr);
623   ierr = PetscOptionsBool("-recursive_split","Split the domain recursively","",ctx.recursive,&ctx.recursive,NULL);CHKERRQ(ierr);
624   ierr = PetscOptionsEnd();CHKERRQ(ierr);
625 
626   /* Choose the limiter from the list of registered limiters */
627   ierr = PetscFunctionListFind(limiters,lname,&ctx.limit);CHKERRQ(ierr);
628   if (!ctx.limit) SETERRQ1(PETSC_COMM_SELF,1,"Limiter '%s' not found",lname);
629 
630   /* Choose the physics from the list of registered models */
631   {
632     PetscErrorCode (*r)(FVCtx*);
633     ierr = PetscFunctionListFind(physics,physname,&r);CHKERRQ(ierr);
634     if (!r) SETERRQ1(PETSC_COMM_SELF,1,"Physics '%s' not found",physname);
635     /* Create the physics, will set the number of fields and their names */
636     ierr = (*r)(&ctx);CHKERRQ(ierr);
637   }
638 
639   /* Create a DMDA to manage the parallel grid */
640   ierr = DMDACreate1d(comm,DM_BOUNDARY_PERIODIC,50,ctx.physics.dof,2,NULL,&da);CHKERRQ(ierr);
641   ierr = DMSetFromOptions(da);CHKERRQ(ierr);
642   ierr = DMSetUp(da);CHKERRQ(ierr);
643   /* Inform the DMDA of the field names provided by the physics. */
644   /* The names will be shown in the title bars when run with -ts_monitor_draw_solution */
645   for (i=0; i<ctx.physics.dof; i++) {
646     ierr = DMDASetFieldName(da,i,ctx.physics.fieldname[i]);CHKERRQ(ierr);
647   }
648   ierr = DMDAGetInfo(da,0, &Mx,0,0, 0,0,0, &dof,0,0,0,0,0);CHKERRQ(ierr);
649   ierr = DMDAGetCorners(da,&xs,0,0,&xm,0,0);CHKERRQ(ierr);
650 
651   /* Set coordinates of cell centers */
652   ierr = DMDASetUniformCoordinates(da,ctx.xmin+0.5*(ctx.xmax-ctx.xmin)/Mx,ctx.xmax+0.5*(ctx.xmax-ctx.xmin)/Mx,0,0,0,0);CHKERRQ(ierr);
653 
654   /* Allocate work space for the Finite Volume solver (so it doesn't have to be reallocated on each function evaluation) */
655   ierr = PetscMalloc4(dof*dof,&ctx.R,dof*dof,&ctx.Rinv,2*dof,&ctx.cjmpLR,1*dof,&ctx.cslope);CHKERRQ(ierr);
656   ierr = PetscMalloc3(2*dof,&ctx.uLR,dof,&ctx.flux,dof,&ctx.speeds);CHKERRQ(ierr);
657 
658   /* Create a vector to store the solution and to save the initial state */
659   ierr = DMCreateGlobalVector(da,&X);CHKERRQ(ierr);
660   ierr = VecDuplicate(X,&X0);CHKERRQ(ierr);
661   ierr = VecDuplicate(X,&R);CHKERRQ(ierr);
662 
663   /*-------------------------------- create index for slow parts and fast parts ----------------------------------------*/
664   ierr = PetscMalloc1(xm*dof,&index_slow);CHKERRQ(ierr);
665   ierr = PetscMalloc1(xm*dof,&index_fast);CHKERRQ(ierr);
666   for (i=xs; i<xs+xm; i++) {
667     if (ctx.xmin+i*(ctx.xmax-ctx.xmin)/(PetscReal)Mx+0.5*(ctx.xmax-ctx.xmin)/(PetscReal)Mx < 0)
668       for (k=0; k<dof; k++) index_slow[islow++] = i*dof+k;
669     else
670       for (k=0; k<dof; k++) index_fast[ifast++] = i*dof+k;
671   }  /* this step need to be changed based on discontinuous point of a */
672   ierr = ISCreateGeneral(PETSC_COMM_WORLD,islow,index_slow,PETSC_COPY_VALUES,&ctx.iss);CHKERRQ(ierr);
673   ierr = ISCreateGeneral(PETSC_COMM_WORLD,ifast,index_fast,PETSC_COPY_VALUES,&ctx.isf);CHKERRQ(ierr);
674 
675   /* Create a time-stepping object */
676   ierr = TSCreate(comm,&ts);CHKERRQ(ierr);
677   ierr = TSSetDM(ts,da);CHKERRQ(ierr);
678   ierr = TSSetRHSFunction(ts,R,FVRHSFunction,&ctx);CHKERRQ(ierr);
679   ierr = TSRHSSplitSetIS(ts,"slow",ctx.iss);CHKERRQ(ierr);
680   ierr = TSRHSSplitSetIS(ts,"fast",ctx.isf);CHKERRQ(ierr);
681   ierr = TSRHSSplitSetRHSFunction(ts,"slow",NULL,FVRHSFunctionslow,&ctx);CHKERRQ(ierr);
682   ierr = TSRHSSplitSetRHSFunction(ts,"fast",NULL,FVRHSFunctionfast,&ctx);CHKERRQ(ierr);
683 
684   if (ctx.recursive) {
685     TS subts;
686     islow = 0;
687     ifast = 0;
688     for (i=xs; i<xs+xm; i++) {
689       PetscReal coord = ctx.xmin+i*(ctx.xmax-ctx.xmin)/(PetscReal)Mx+0.5*(ctx.xmax-ctx.xmin)/(PetscReal)Mx;
690       if (coord >= 0 && coord < ctx.xmin+(ctx.xmax-ctx.xmin)*3/4.)
691         for (k=0; k<dof; k++) index_slow[islow++] = i*dof+k;
692       if (coord >= ctx.xmin+(ctx.xmax-ctx.xmin)*3/4.)
693         for (k=0; k<dof; k++) index_fast[ifast++] = i*dof+k;
694     }  /* this step need to be changed based on discontinuous point of a */
695     ierr = ISCreateGeneral(PETSC_COMM_WORLD,islow,index_slow,PETSC_COPY_VALUES,&ctx.iss2);CHKERRQ(ierr);
696     ierr = ISCreateGeneral(PETSC_COMM_WORLD,ifast,index_fast,PETSC_COPY_VALUES,&ctx.isf2);CHKERRQ(ierr);
697 
698     ierr = TSRHSSplitGetSubTS(ts,"fast",&subts);CHKERRQ(ierr);
699     ierr = TSRHSSplitSetIS(subts,"slow",ctx.iss2);CHKERRQ(ierr);
700     ierr = TSRHSSplitSetIS(subts,"fast",ctx.isf2);CHKERRQ(ierr);
701     ierr = TSRHSSplitSetRHSFunction(subts,"slow",NULL,FVRHSFunctionslow2,&ctx);CHKERRQ(ierr);
702     ierr = TSRHSSplitSetRHSFunction(subts,"fast",NULL,FVRHSFunctionfast2,&ctx);CHKERRQ(ierr);
703   }
704 
705   /*ierr = TSSetType(ts,TSSSP);CHKERRQ(ierr);*/
706   ierr = TSSetType(ts,TSMPRK);CHKERRQ(ierr);
707   ierr = TSSetMaxTime(ts,10);CHKERRQ(ierr);
708   ierr = TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);CHKERRQ(ierr);
709 
710   /* Compute initial conditions and starting time step */
711   ierr = FVSample(&ctx,da,0,X0);CHKERRQ(ierr);
712   ierr = FVRHSFunction(ts,0,X0,X,(void*)&ctx);CHKERRQ(ierr); /* Initial function evaluation, only used to determine max speed */
713   ierr = VecCopy(X0,X);CHKERRQ(ierr);                        /* The function value was not used so we set X=X0 again */
714   ierr = TSSetTimeStep(ts,ctx.cfl/ctx.cfl_idt);CHKERRQ(ierr);
715   ierr = TSSetFromOptions(ts);CHKERRQ(ierr); /* Take runtime options */
716   ierr = SolutionStatsView(da,X,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
717   {
718     PetscInt    steps;
719     PetscScalar mass_initial,mass_final,mass_difference;
720 
721     ierr = TSSolve(ts,X);CHKERRQ(ierr);
722     ierr = TSGetSolveTime(ts,&ptime);CHKERRQ(ierr);
723     ierr = TSGetStepNumber(ts,&steps);CHKERRQ(ierr);
724     ierr = PetscPrintf(comm,"Final time %g, steps %D\n",(double)ptime,steps);CHKERRQ(ierr);
725     /* calculate the total mass at initial time and final time */
726     mass_initial = 0.0;
727     mass_final   = 0.0;
728     ierr = VecSum(X0,&mass_initial);CHKERRQ(ierr);
729     ierr = VecSum(X,&mass_final);CHKERRQ(ierr);
730     mass_difference = (ctx.xmax-ctx.xmin)/(PetscScalar)Mx*(mass_final - mass_initial);
731     ierr = PetscPrintf(comm,"Mass difference %g\n",(double)mass_difference);CHKERRQ(ierr);
732     if (ctx.simulation) {
733       PetscViewer  fd;
734       char         filename[PETSC_MAX_PATH_LEN] = "binaryoutput";
735       Vec          XR;
736       PetscReal    nrm1,nrmsup;
737       PetscBool    flg;
738       ierr = PetscOptionsGetString(NULL,NULL,"-f",filename,sizeof(filename),&flg);CHKERRQ(ierr);
739       if (!flg) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_USER,"Must indicate binary file with the -f option");
740       ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_READ,&fd);CHKERRQ(ierr);
741       ierr = VecDuplicate(X0,&XR);CHKERRQ(ierr);
742       ierr = VecLoad(XR,fd);CHKERRQ(ierr);
743       ierr = PetscViewerDestroy(&fd);CHKERRQ(ierr);
744       ierr = VecAYPX(XR,-1,X);CHKERRQ(ierr);
745       ierr = VecNorm(XR,NORM_1,&nrm1);CHKERRQ(ierr);
746       ierr = VecNorm(XR,NORM_INFINITY,&nrmsup);CHKERRQ(ierr);
747       nrm1 /= Mx;
748       ierr = PetscPrintf(comm,"Error ||x-x_e||_1 %g  ||x-x_e||_sup %g\n",(double)nrm1,(double)nrmsup);CHKERRQ(ierr);
749       ierr = VecDestroy(&XR);CHKERRQ(ierr);
750     }
751   }
752 
753   ierr = SolutionStatsView(da,X,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
754   if (draw & 0x1) {ierr = VecView(X0,PETSC_VIEWER_DRAW_WORLD);CHKERRQ(ierr);}
755   if (draw & 0x2) {ierr = VecView(X,PETSC_VIEWER_DRAW_WORLD);CHKERRQ(ierr);}
756   if (draw & 0x4) {
757     Vec Y;
758     ierr = VecDuplicate(X,&Y);CHKERRQ(ierr);
759     ierr = FVSample(&ctx,da,ptime,Y);CHKERRQ(ierr);
760     ierr = VecAYPX(Y,-1,X);CHKERRQ(ierr);
761     ierr = VecView(Y,PETSC_VIEWER_DRAW_WORLD);CHKERRQ(ierr);
762     ierr = VecDestroy(&Y);CHKERRQ(ierr);
763   }
764 
765   if (view_final) {
766     PetscViewer viewer;
767     ierr = PetscViewerASCIIOpen(PETSC_COMM_WORLD,final_fname,&viewer);CHKERRQ(ierr);
768     ierr = PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_MATLAB);CHKERRQ(ierr);
769     ierr = VecView(X,viewer);CHKERRQ(ierr);
770     ierr = PetscViewerPopFormat(viewer);CHKERRQ(ierr);
771     ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);
772   }
773 
774   /* Clean up */
775   ierr = (*ctx.physics.destroy)(ctx.physics.user);CHKERRQ(ierr);
776   for (i=0; i<ctx.physics.dof; i++) {ierr = PetscFree(ctx.physics.fieldname[i]);CHKERRQ(ierr);}
777   ierr = PetscFree4(ctx.R,ctx.Rinv,ctx.cjmpLR,ctx.cslope);CHKERRQ(ierr);
778   ierr = PetscFree3(ctx.uLR,ctx.flux,ctx.speeds);CHKERRQ(ierr);
779   ierr = ISDestroy(&ctx.iss);CHKERRQ(ierr);
780   ierr = ISDestroy(&ctx.isf);CHKERRQ(ierr);
781   ierr = ISDestroy(&ctx.iss2);CHKERRQ(ierr);
782   ierr = ISDestroy(&ctx.isf2);CHKERRQ(ierr);
783   ierr = VecDestroy(&X);CHKERRQ(ierr);
784   ierr = VecDestroy(&X0);CHKERRQ(ierr);
785   ierr = VecDestroy(&R);CHKERRQ(ierr);
786   ierr = DMDestroy(&da);CHKERRQ(ierr);
787   ierr = TSDestroy(&ts);CHKERRQ(ierr);
788   ierr = PetscFree(index_slow);CHKERRQ(ierr);
789   ierr = PetscFree(index_fast);CHKERRQ(ierr);
790   ierr = PetscFunctionListDestroy(&limiters);CHKERRQ(ierr);
791   ierr = PetscFunctionListDestroy(&physics);CHKERRQ(ierr);
792   ierr = PetscFinalize();
793   return ierr;
794 }
795 
796 /*TEST
797     build:
798       requires: !complex
799       depends: finitevolume1d.c
800 
801     test:
802       args: -da_grid_x 60 -initial 1 -xmin -1 -xmax 1 -limit mc -physics_advect_a1 1 -physics_advect_a2 2 -ts_dt 0.025 -ts_max_steps 24 -ts_type rk -ts_rk_type 2a -ts_rk_dtratio 2 -ts_rk_multirate -ts_use_splitrhsfunction 0
803 
804     test:
805       suffix: 2
806       args: -da_grid_x 60 -initial 1 -xmin -1 -xmax 1 -limit mc -physics_advect_a1 1 -physics_advect_a2 2 -ts_dt 0.025 -ts_max_steps 24 -ts_type rk -ts_rk_type 2a -ts_rk_dtratio 2 -ts_rk_multirate -ts_use_splitrhsfunction 1
807       output_file: output/ex5_1.out
808 
809     test:
810       suffix: 3
811       args: -da_grid_x 60 -initial 1 -xmin -1 -xmax 1 -limit mc -physics_advect_a1 1 -physics_advect_a2 2 -ts_dt 0.025 -ts_max_steps 24 -ts_type mprk -ts_mprk_type 2a22 -ts_use_splitrhsfunction 0
812 
813     test:
814       suffix: 4
815       args: -da_grid_x 60 -initial 1 -xmin -1 -xmax 1 -limit mc -physics_advect_a1 1 -physics_advect_a2 2 -ts_dt 0.025 -ts_max_steps 24 -ts_type mprk -ts_mprk_type 2a22 -ts_use_splitrhsfunction 1
816       output_file: output/ex5_3.out
817 
818     test:
819       suffix: 5
820       args: -da_grid_x 60 -initial 1 -xmin -1 -xmax 1 -limit mc -physics_advect_a1 1 -physics_advect_a2 2 -ts_dt 0.025 -ts_max_steps 24 -ts_type rk -ts_rk_type 2a -ts_rk_dtratio 2 -ts_rk_multirate -ts_use_splitrhsfunction 0 -recursive_split
821 
822     test:
823       suffix: 6
824       args: -da_grid_x 60 -initial 1 -xmin -1 -xmax 1 -limit mc -physics_advect_a1 1 -physics_advect_a2 2 -ts_dt 0.025 -ts_max_steps 24 -ts_type rk -ts_rk_type 2a -ts_rk_dtratio 2 -ts_rk_multirate -ts_use_splitrhsfunction 1 -recursive_split
825 TEST*/
826