Actual source code: ex7.c
1: static const char help[] = "1D periodic Finite Volume solver by a particular slope limiter with semidiscrete time stepping.\n"
2: " advection - Constant coefficient scalar advection\n"
3: " u_t + (a*u)_x = 0\n"
4: " for this toy problem, we choose different meshsizes for different sub-domains, say\n"
5: " hxs = (xmax - xmin)/2.0*(hratio+1.0)/Mx, \n"
6: " hxf = (xmax - xmin)/2.0*(1.0+1.0/hratio)/Mx, \n"
7: " with x belongs to (xmin,xmax), the number of total mesh points is Mx and the ratio between the meshsize of corse\n\n"
8: " grids and fine grids is hratio.\n"
9: " exact - Exact Riemann solver which usually needs to perform a Newton iteration to connect\n"
10: " the states across shocks and rarefactions\n"
11: " simulation - use reference solution which is generated by smaller time step size to be true solution,\n"
12: " also the reference solution should be generated by user and stored in a binary file.\n"
13: " characteristic - Limit the characteristic variables, this is usually preferred (default)\n"
14: "Several initial conditions can be chosen with -initial N\n\n"
15: "The problem size should be set with -da_grid_x M\n\n"
16: "This script choose the slope limiter by biased second-order upwind procedure which is proposed by Van Leer in 1994\n"
17: " u(x_(k+1/2),t) = u(x_k,t) + phi(x_(k+1/2),t)*(u(x_k,t)-u(x_(k-1),t)) \n"
18: " limiter phi(x_(k+1/2),t) = max(0,min(r(k+1/2),min(2,gamma(k+1/2)*r(k+1/2)+alpha(k+1/2)))) \n"
19: " r(k+1/2) = (u(x_(k+1))-u(x_k))/(u(x_k)-u(x_(k-1))) \n"
20: " alpha(k+1/2) = (h_k*h_(k+1))/(h_(k-1)+h_k)/(h_(k-1)+h_k+h_(k+1)) \n"
21: " gamma(k+1/2) = h_k*(h_(k-1)+h_k)/(h_k+h_(k+1))/(h_(k-1)+h_k+h_(k+1)) \n";
23: #include <petscts.h>
24: #include <petscdm.h>
25: #include <petscdmda.h>
26: #include <petscdraw.h>
27: #include <petscmath.h>
29: PETSC_STATIC_INLINE PetscReal RangeMod(PetscReal a,PetscReal xmin,PetscReal xmax) { PetscReal range = xmax-xmin; return xmin +PetscFmodReal(range+PetscFmodReal(a,range),range); }
31: /* --------------------------------- Finite Volume data structures ----------------------------------- */
33: typedef enum {FVBC_PERIODIC, FVBC_OUTFLOW} FVBCType;
34: static const char *FVBCTypes[] = {"PERIODIC","OUTFLOW","FVBCType","FVBC_",0};
36: typedef struct {
37: PetscErrorCode (*sample)(void*,PetscInt,FVBCType,PetscReal,PetscReal,PetscReal,PetscReal,PetscReal*);
38: PetscErrorCode (*flux)(void*,const PetscScalar*,PetscScalar*,PetscReal*);
39: PetscErrorCode (*destroy)(void*);
40: void *user;
41: PetscInt dof;
42: char *fieldname[16];
43: } PhysicsCtx;
45: typedef struct {
46: PhysicsCtx physics;
47: MPI_Comm comm;
48: char prefix[256];
50: /* Local work arrays */
51: PetscScalar *flux; /* Flux across interface */
52: PetscReal *speeds; /* Speeds of each wave */
53: PetscScalar *u; /* value at face */
55: PetscReal cfl_idt; /* Max allowable value of 1/Delta t */
56: PetscReal cfl;
57: PetscReal xmin,xmax;
58: PetscInt initial;
59: PetscBool exact;
60: PetscBool simulation;
61: FVBCType bctype;
62: PetscInt hratio; /* hratio = hslow/hfast */
63: IS isf,iss;
64: PetscInt sf,fs; /* slow-fast and fast-slow interfaces */
65: } FVCtx;
67: /* --------------------------------- Physics ----------------------------------- */
68: static PetscErrorCode PhysicsDestroy_SimpleFree(void *vctx)
69: {
73: PetscFree(vctx);
74: return(0);
75: }
77: /* --------------------------------- Advection ----------------------------------- */
78: typedef struct {
79: PetscReal a; /* advective velocity */
80: } AdvectCtx;
82: static PetscErrorCode PhysicsFlux_Advect(void *vctx,const PetscScalar *u,PetscScalar *flux,PetscReal *maxspeed)
83: {
84: AdvectCtx *ctx = (AdvectCtx*)vctx;
85: PetscReal speed;
88: speed = ctx->a;
89: flux[0] = speed*u[0];
90: *maxspeed = speed;
91: return(0);
92: }
94: static PetscErrorCode PhysicsSample_Advect(void *vctx,PetscInt initial,FVBCType bctype,PetscReal xmin,PetscReal xmax,PetscReal t,PetscReal x,PetscReal *u)
95: {
96: AdvectCtx *ctx = (AdvectCtx*)vctx;
97: PetscReal a = ctx->a,x0;
100: switch (bctype) {
101: case FVBC_OUTFLOW: x0 = x-a*t; break;
102: case FVBC_PERIODIC: x0 = RangeMod(x-a*t,xmin,xmax); break;
103: default: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"unknown BCType");
104: }
105: switch (initial) {
106: case 0: u[0] = (x0 < 0) ? 1 : -1; break;
107: case 1: u[0] = (x0 < 0) ? -1 : 1; break;
108: case 2: u[0] = (0 < x0 && x0 < 1) ? 1 : 0; break;
109: case 3: u[0] = PetscSinReal(2*PETSC_PI*x0); break;
110: case 4: u[0] = PetscAbs(x0); break;
111: case 5: u[0] = (x0 < 0 || x0 > 0.5) ? 0 : PetscSqr(PetscSinReal(2*PETSC_PI*x0)); break;
112: case 6: u[0] = (x0 < 0) ? 0 : ((x0 < 1) ? x0 : ((x0 < 2) ? 2-x0 : 0)); break;
113: case 7: u[0] = PetscPowReal(PetscSinReal(PETSC_PI*x0),10.0);break;
114: default: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"unknown initial condition");
115: }
116: return(0);
117: }
119: static PetscErrorCode PhysicsCreate_Advect(FVCtx *ctx)
120: {
122: AdvectCtx *user;
125: PetscNew(&user);
126: ctx->physics.sample = PhysicsSample_Advect;
127: ctx->physics.flux = PhysicsFlux_Advect;
128: ctx->physics.destroy = PhysicsDestroy_SimpleFree;
129: ctx->physics.user = user;
130: ctx->physics.dof = 1;
131: PetscStrallocpy("u",&ctx->physics.fieldname[0]);
132: user->a = 1;
133: PetscOptionsBegin(ctx->comm,ctx->prefix,"Options for advection","");
134: {
135: PetscOptionsReal("-physics_advect_a","Speed","",user->a,&user->a,NULL);
136: }
137: PetscOptionsEnd();
138: return(0);
139: }
141: /* --------------------------------- Finite Volume Solver ----------------------------------- */
143: static PetscErrorCode FVRHSFunction(TS ts,PetscReal time,Vec X,Vec F,void *vctx)
144: {
145: FVCtx *ctx = (FVCtx*)vctx;
147: PetscInt i,j,Mx,dof,xs,xm,sf = ctx->sf,fs = ctx->fs;
148: PetscReal hf,hs,cfl_idt = 0;
149: PetscScalar *x,*f,*r,*min,*alpha,*gamma;
150: Vec Xloc;
151: DM da;
154: TSGetDM(ts,&da);
155: DMGetLocalVector(da,&Xloc); /* Xloc contains ghost points */
156: DMDAGetInfo(da,0,&Mx,0,0,0,0,0,&dof,0,0,0,0,0); /* Mx is the number of center points */
157: hs = (ctx->xmax-ctx->xmin)/2.0*(ctx->hratio+1.0)/Mx;
158: hf = (ctx->xmax-ctx->xmin)/2.0*(1.0+1.0/ctx->hratio)/Mx;
159: DMGlobalToLocalBegin(da,X,INSERT_VALUES,Xloc); /* X is solution vector which does not contain ghost points */
160: DMGlobalToLocalEnd(da,X,INSERT_VALUES,Xloc);
161: VecZeroEntries(F); /* F is the right hand side function corresponds to center points */
162: DMDAVecGetArray(da,Xloc,&x);
163: DMDAVecGetArray(da,F,&f);
164: DMDAGetCorners(da,&xs,0,0,&xm,0,0);
165: PetscMalloc4(dof,&r,dof,&min,dof,&alpha,dof,&gamma);
167: if (ctx->bctype == FVBC_OUTFLOW) {
168: for (i=xs-2; i<0; i++) {
169: for (j=0; j<dof; j++) x[i*dof+j] = x[j];
170: }
171: for (i=Mx; i<xs+xm+2; i++) {
172: for (j=0; j<dof; j++) x[i*dof+j] = x[(xs+xm-1)*dof+j];
173: }
174: }
176: for (i=xs; i<xs+xm+1; i++) {
177: PetscReal maxspeed;
178: PetscScalar *u;
179: if (i < sf || i > fs+1) {
180: u = &ctx->u[0];
181: alpha[0] = 1.0/6.0;
182: gamma[0] = 1.0/3.0;
183: for (j=0; j<dof; j++) {
184: r[j] = (x[i*dof+j]-x[(i-1)*dof+j])/(x[(i-1)*dof+j]-x[(i-2)*dof+j]);
185: min[j] = PetscMin(r[j],2.0);
186: u[j] = x[(i-1)*dof+j]+PetscMax(0,PetscMin(min[j],alpha[0]+gamma[0]*r[j]))*(x[(i-1)*dof+j]-x[(i-2)*dof+j]);
187: }
188: (*ctx->physics.flux)(ctx->physics.user,u,ctx->flux,&maxspeed);
189: cfl_idt = PetscMax(cfl_idt,PetscAbsScalar(maxspeed/hs));
190: if (i > xs) {
191: for (j=0; j<dof; j++) f[(i-1)*dof+j] -= ctx->flux[j]/hs;
192: }
193: if (i < xs+xm) {
194: for (j=0; j<dof; j++) f[i*dof+j] += ctx->flux[j]/hs;
195: }
196: } else if (i == sf) {
197: u = &ctx->u[0];
198: alpha[0] = hs*hf/(hs+hs)/(hs+hs+hf);
199: gamma[0] = hs*(hs+hs)/(hs+hf)/(hs+hs+hf);
200: for (j=0; j<dof; j++) {
201: r[j] = (x[i*dof+j]-x[(i-1)*dof+j])/(x[(i-1)*dof+j]-x[(i-2)*dof+j]);
202: min[j] = PetscMin(r[j],2.0);
203: u[j] = x[(i-1)*dof+j]+PetscMax(0,PetscMin(min[j],alpha[0]+gamma[0]*r[j]))*(x[(i-1)*dof+j]-x[(i-2)*dof+j]);
204: }
205: (*ctx->physics.flux)(ctx->physics.user,u,ctx->flux,&maxspeed);
206: if (i > xs) {
207: for (j=0; j<dof; j++) f[(i-1)*dof+j] -= ctx->flux[j]/hs;
208: }
209: if (i < xs+xm) {
210: for (j=0; j<dof; j++) f[i*dof+j] += ctx->flux[j]/hf;
211: }
212: } else if (i == sf+1) {
213: u = &ctx->u[0];
214: alpha[0] = hf*hf/(hs+hf)/(hs+hf+hf);
215: gamma[0] = hf*(hs+hf)/(hf+hf)/(hs+hf+hf);
216: for (j=0; j<dof; j++) {
217: r[j] = (x[i*dof+j] - x[(i-1)*dof+j])/(x[(i-1)*dof+j]-x[(i-2)*dof+j]);
218: min[j] = PetscMin(r[j],2.0);
219: u[j] = x[(i-1)*dof+j]+PetscMax(0,PetscMin(min[j],alpha[0]+gamma[0]*r[j]))*(x[(i-1)*dof+j]-x[(i-2)*dof+j]);
220: }
221: (*ctx->physics.flux)(ctx->physics.user,u,ctx->flux,&maxspeed);
222: if (i > xs) {
223: for (j=0; j<dof; j++) f[(i-1)*dof+j] -= ctx->flux[j]/hf;
224: }
225: if (i < xs+xm) {
226: for (j=0; j<dof; j++) f[i*dof+j] += ctx->flux[j]/hf;
227: }
228: } else if (i > sf+1 && i < fs) {
229: u = &ctx->u[0];
230: alpha[0] = 1.0/6.0;
231: gamma[0] = 1.0/3.0;
232: for (j=0; j<dof; j++) {
233: r[j] = (x[i*dof+j] - x[(i-1)*dof+j])/(x[(i-1)*dof+j]-x[(i-2)*dof+j]);
234: min[j] = PetscMin(r[j],2.0);
235: u[j] = x[(i-1)*dof+j]+PetscMax(0,PetscMin(min[j],alpha[0]+gamma[0]*r[j]))*(x[(i-1)*dof+j]-x[(i-2)*dof+j]);
236: }
237: (*ctx->physics.flux)(ctx->physics.user,u,ctx->flux,&maxspeed);
238: if (i > xs) {
239: for (j=0; j<dof; j++) f[(i-1)*dof+j] -= ctx->flux[j]/hf;
240: }
241: if (i < xs+xm) {
242: for (j=0; j<dof; j++) f[i*dof+j] += ctx->flux[j]/hf;
243: }
244: } else if (i == fs) {
245: u = &ctx->u[0];
246: alpha[0] = hf*hs/(hf+hf)/(hf+hf+hs);
247: gamma[0] = hf*(hf+hf)/(hf+hs)/(hf+hf+hs);
248: for (j=0; j<dof; j++) {
249: r[j] = (x[i*dof+j] - x[(i-1)*dof+j])/(x[(i-1)*dof+j]-x[(i-2)*dof+j]);
250: min[j] = PetscMin(r[j],2.0);
251: u[j] = x[(i-1)*dof+j]+PetscMax(0,PetscMin(min[j],alpha[0]+gamma[0]*r[j]))*(x[(i-1)*dof+j]-x[(i-2)*dof+j]);
252: }
253: (*ctx->physics.flux)(ctx->physics.user,u,ctx->flux,&maxspeed);
254: if (i > xs) {
255: for (j=0; j<dof; j++) f[(i-1)*dof+j] -= ctx->flux[j]/hf;
256: }
257: if (i < xs+xm) {
258: for (j=0; j<dof; j++) f[i*dof+j] += ctx->flux[j]/hs;
259: }
260: } else if (i == fs+1) {
261: u = &ctx->u[0];
262: alpha[0] = hs*hs/(hf+hs)/(hf+hs+hs);
263: gamma[0] = hs*(hf+hs)/(hs+hs)/(hf+hs+hs);
264: for (j=0; j<dof; j++) {
265: r[j] = (x[i*dof+j] - x[(i-1)*dof+j])/(x[(i-1)*dof+j]-x[(i-2)*dof+j]);
266: min[j] = PetscMin(r[j],2.0);
267: u[j] = x[(i-1)*dof+j]+PetscMax(0,PetscMin(min[j],alpha[0]+gamma[0]*r[j]))*(x[(i-1)*dof+j]-x[(i-2)*dof+j]);
268: }
269: (*ctx->physics.flux)(ctx->physics.user,u,ctx->flux,&maxspeed);
270: if (i > xs) {
271: for (j=0; j<dof; j++) f[(i-1)*dof+j] -= ctx->flux[j]/hs;
272: }
273: if (i < xs+xm) {
274: for (j=0; j<dof; j++) f[i*dof+j] += ctx->flux[j]/hs;
275: }
276: }
277: }
278: DMDAVecRestoreArray(da,Xloc,&x);
279: DMDAVecRestoreArray(da,F,&f);
280: DMRestoreLocalVector(da,&Xloc);
281: MPI_Allreduce(&cfl_idt,&ctx->cfl_idt,1,MPIU_SCALAR,MPIU_MAX,PetscObjectComm((PetscObject)da));
282: if (0) {
283: /* We need a way to inform the TS of a CFL constraint, this is a debugging fragment */
284: PetscReal dt,tnow;
285: TSGetTimeStep(ts,&dt);
286: TSGetTime(ts,&tnow);
287: if (dt > 0.5/ctx->cfl_idt) {
288: PetscPrintf(ctx->comm,"Stability constraint exceeded at t=%g, dt %g > %g\n",(double)tnow,(double)dt,(double)(0.5/ctx->cfl_idt));
289: }
290: }
291: PetscFree4(r,min,alpha,gamma);
292: return(0);
293: }
295: static PetscErrorCode FVRHSFunctionslow(TS ts,PetscReal time,Vec X,Vec F,void *vctx)
296: {
297: FVCtx *ctx = (FVCtx*)vctx;
298: PetscErrorCode ierr;
299: PetscInt i,j,Mx,dof,xs,xm,islow = 0,sf = ctx->sf,fs = ctx->fs;
300: PetscReal hf,hs;
301: PetscScalar *x,*f,*r,*min,*alpha,*gamma;
302: Vec Xloc;
303: DM da;
306: TSGetDM(ts,&da);
307: DMGetLocalVector(da,&Xloc); /* Xloc contains ghost points */
308: DMDAGetInfo(da,0,&Mx,0,0,0,0,0,&dof,0,0,0,0,0); /* Mx is the number of center points */
309: hs = (ctx->xmax-ctx->xmin)/2.0*(ctx->hratio+1.0)/Mx;
310: hf = (ctx->xmax-ctx->xmin)/2.0*(1.0+1.0/ctx->hratio)/Mx;
311: DMGlobalToLocalBegin(da,X,INSERT_VALUES,Xloc); /* X is solution vector which does not contain ghost points */
312: DMGlobalToLocalEnd (da,X,INSERT_VALUES,Xloc);
313: VecZeroEntries(F); /* F is the right hand side function corresponds to center points */
314: DMDAVecGetArray(da,Xloc,&x);
315: VecGetArray(F,&f);
316: DMDAGetCorners(da,&xs,0,0,&xm,0,0);
317: PetscMalloc4(dof,&r,dof,&min,dof,&alpha,dof,&gamma);
319: if (ctx->bctype == FVBC_OUTFLOW) {
320: for (i=xs-2; i<0; i++) {
321: for (j=0; j<dof; j++) x[i*dof+j] = x[j];
322: }
323: for (i=Mx; i<xs+xm+2; i++) {
324: for (j=0; j<dof; j++) x[i*dof+j] = x[(xs+xm-1)*dof+j];
325: }
326: }
328: for (i=xs; i<xs+xm+1; i++) {
329: PetscReal maxspeed;
330: PetscScalar *u;
331: if (i < sf) {
332: u = &ctx->u[0];
333: alpha[0] = 1.0/6.0;
334: gamma[0] = 1.0/3.0;
335: for (j=0; j<dof; j++) {
336: r[j] = (x[i*dof+j] - x[(i-1)*dof+j])/(x[(i-1)*dof+j]-x[(i-2)*dof+j]);
337: min[j] = PetscMin(r[j],2.0);
338: u[j] = x[(i-1)*dof+j]+PetscMax(0,PetscMin(min[j],alpha[0]+gamma[0]*r[j]))*(x[(i-1)*dof+j]-x[(i-2)*dof+j]);
339: }
340: (*ctx->physics.flux)(ctx->physics.user,u,ctx->flux,&maxspeed);
341: if (i > xs) {
342: for (j=0; j<dof; j++) f[(islow-1)*dof+j] -= ctx->flux[j]/hs;
343: }
344: if (i < xs+xm) {
345: for (j=0; j<dof; j++) f[islow*dof+j] += ctx->flux[j]/hs;
346: islow++;
347: }
348: } else if (i == sf) {
349: u = &ctx->u[0];
350: alpha[0] = hs*hf/(hs+hs)/(hs+hs+hf);
351: gamma[0] = hs*(hs+hs)/(hs+hf)/(hs+hs+hf);
352: for (j=0; j<dof; j++) {
353: r[j] = (x[i*dof+j] - x[(i-1)*dof+j])/(x[(i-1)*dof+j]-x[(i-2)*dof+j]);
354: min[j] = PetscMin(r[j],2.0);
355: u[j] = x[(i-1)*dof+j]+PetscMax(0,PetscMin(min[j],alpha[0]+gamma[0]*r[j]))*(x[(i-1)*dof+j]-x[(i-2)*dof+j]);
356: }
357: (*ctx->physics.flux)(ctx->physics.user,u,ctx->flux,&maxspeed);
358: if (i > xs) {
359: for (j=0; j<dof; j++) f[(islow-1)*dof+j] -= ctx->flux[j]/hs;
360: }
361: } else if (i == fs) {
362: u = &ctx->u[0];
363: alpha[0] = hf*hs/(hf+hf)/(hf+hf+hs);
364: gamma[0] = hf*(hf+hf)/(hf+hs)/(hf+hf+hs);
365: for (j=0; j<dof; j++) {
366: r[j] = (x[i*dof+j] - x[(i-1)*dof+j])/(x[(i-1)*dof+j]-x[(i-2)*dof+j]);
367: min[j] = PetscMin(r[j],2.0);
368: u[j] = x[(i-1)*dof+j]+PetscMax(0,PetscMin(min[j],alpha[0]+gamma[0]*r[j]))*(x[(i-1)*dof+j]-x[(i-2)*dof+j]);
369: }
370: (*ctx->physics.flux)(ctx->physics.user,u,ctx->flux,&maxspeed);
371: if (i < xs+xm) {
372: for (j=0; j<dof; j++) f[islow*dof+j] += ctx->flux[j]/hs;
373: islow++;
374: }
375: } else if (i == fs+1) {
376: u = &ctx->u[0];
377: alpha[0] = hs*hs/(hf+hs)/(hf+hs+hs);
378: gamma[0] = hs*(hf+hs)/(hs+hs)/(hf+hs+hs);
379: for (j=0; j<dof; j++) {
380: r[j] = (x[i*dof+j] - x[(i-1)*dof+j])/(x[(i-1)*dof+j]-x[(i-2)*dof+j]);
381: min[j] = PetscMin(r[j],2.0);
382: u[j] = x[(i-1)*dof+j]+PetscMax(0,PetscMin(min[j],alpha[0]+gamma[0]*r[j]))*(x[(i-1)*dof+j]-x[(i-2)*dof+j]);
383: }
384: (*ctx->physics.flux)(ctx->physics.user,u,ctx->flux,&maxspeed);
385: if (i > xs) {
386: for (j=0; j<dof; j++) f[(islow-1)*dof+j] -= ctx->flux[j]/hs;
387: }
388: if (i < xs+xm) {
389: for (j=0; j<dof; j++) f[islow*dof+j] += ctx->flux[j]/hs;
390: islow++;
391: }
392: } else if (i > fs+1) {
393: u = &ctx->u[0];
394: alpha[0] = 1.0/6.0;
395: gamma[0] = 1.0/3.0;
396: for (j=0; j<dof; j++) {
397: r[j] = (x[i*dof+j] - x[(i-1)*dof+j])/(x[(i-1)*dof+j]-x[(i-2)*dof+j]);
398: min[j] = PetscMin(r[j],2.0);
399: u[j] = x[(i-1)*dof+j]+PetscMax(0,PetscMin(min[j],alpha[0]+gamma[0]*r[j]))*(x[(i-1)*dof+j]-x[(i-2)*dof+j]);
400: }
401: (*ctx->physics.flux)(ctx->physics.user,u,ctx->flux,&maxspeed);
402: if (i > xs) {
403: for (j=0; j<dof; j++) f[(islow-1)*dof+j] -= ctx->flux[j]/hs;
404: }
405: if (i < xs+xm) {
406: for (j=0; j<dof; j++) f[islow*dof+j] += ctx->flux[j]/hs;
407: islow++;
408: }
409: }
410: }
411: DMDAVecRestoreArray(da,Xloc,&x);
412: VecRestoreArray(F,&f);
413: DMRestoreLocalVector(da,&Xloc);
414: PetscFree4(r,min,alpha,gamma);
415: return(0);
416: }
418: static PetscErrorCode FVRHSFunctionfast(TS ts,PetscReal time,Vec X,Vec F,void *vctx)
419: {
420: FVCtx *ctx = (FVCtx*)vctx;
422: PetscInt i,j,Mx,dof,xs,xm,ifast = 0,sf = ctx->sf,fs = ctx->fs;
423: PetscReal hf,hs;
424: PetscScalar *x,*f,*r,*min,*alpha,*gamma;
425: Vec Xloc;
426: DM da;
429: TSGetDM(ts,&da);
430: DMGetLocalVector(da,&Xloc); /* Xloc contains ghost points */
431: DMDAGetInfo(da,0,&Mx,0,0,0,0,0,&dof,0,0,0,0,0); /* Mx is the number of center points */
432: hs = (ctx->xmax-ctx->xmin)/2.0*(ctx->hratio+1.0)/Mx;
433: hf = (ctx->xmax-ctx->xmin)/2.0*(1.0+1.0/ctx->hratio)/Mx;
434: DMGlobalToLocalBegin(da,X,INSERT_VALUES,Xloc); /* X is solution vector which does not contain ghost points */
435: DMGlobalToLocalEnd(da,X,INSERT_VALUES,Xloc);
436: VecZeroEntries(F); /* F is the right hand side function corresponds to center points */
437: DMDAVecGetArray(da,Xloc,&x);
438: VecGetArray(F,&f);
439: DMDAGetCorners(da,&xs,0,0,&xm,0,0);
440: PetscMalloc4(dof,&r,dof,&min,dof,&alpha,dof,&gamma);
442: if (ctx->bctype == FVBC_OUTFLOW) {
443: for (i=xs-2; i<0; i++) {
444: for (j=0; j<dof; j++) x[i*dof+j] = x[j];
445: }
446: for (i=Mx; i<xs+xm+2; i++) {
447: for (j=0; j<dof; j++) x[i*dof+j] = x[(xs+xm-1)*dof+j];
448: }
449: }
451: for (i=xs; i<xs+xm+1; i++) {
452: PetscReal maxspeed;
453: PetscScalar *u;
454: if (i == sf) {
455: u = &ctx->u[0];
456: alpha[0] = hs*hf/(hs+hs)/(hs+hs+hf);
457: gamma[0] = hs*(hs+hs)/(hs+hf)/(hs+hs+hf);
458: for (j=0; j<dof; j++) {
459: r[j] = (x[i*dof+j]-x[(i-1)*dof+j])/(x[(i-1)*dof+j]-x[(i-2)*dof+j]);
460: min[j] = PetscMin(r[j],2.0);
461: u[j] = x[(i-1)*dof+j]+PetscMax(0,PetscMin(min[j],alpha[0]+gamma[0]*r[j]))*(x[(i-1)*dof+j]-x[(i-2)*dof+j]);
462: }
463: (*ctx->physics.flux)(ctx->physics.user,u,ctx->flux,&maxspeed);
464: if (i < xs+xm) {
465: for (j=0; j<dof; j++) f[ifast*dof+j] += ctx->flux[j]/hf;
466: ifast++;
467: }
468: } else if (i == sf+1) {
469: u = &ctx->u[0];
470: alpha[0] = hf*hf/(hs+hf)/(hs+hf+hf);
471: gamma[0] = hf*(hs+hf)/(hf+hf)/(hs+hf+hf);
472: for (j=0; j<dof; j++) {
473: r[j] = (x[i*dof+j]-x[(i-1)*dof+j])/(x[(i-1)*dof+j]-x[(i-2)*dof+j]);
474: min[j] = PetscMin(r[j],2.0);
475: u[j] = x[(i-1)*dof+j]+PetscMax(0,PetscMin(min[j],alpha[0]+gamma[0]*r[j]))*(x[(i-1)*dof+j]-x[(i-2)*dof+j]);
476: }
477: (*ctx->physics.flux)(ctx->physics.user,u,ctx->flux,&maxspeed);
478: if (i > xs) {
479: for (j=0; j<dof; j++) f[(ifast-1)*dof+j] -= ctx->flux[j]/hf;
480: }
481: if (i < xs+xm) {
482: for (j=0; j<dof; j++) f[ifast*dof+j] += ctx->flux[j]/hf;
483: ifast++;
484: }
485: } else if (i > sf+1 && i < fs) {
486: u = &ctx->u[0];
487: alpha[0] = 1.0/6.0;
488: gamma[0] = 1.0/3.0;
489: for (j=0; j<dof; j++) {
490: r[j] = (x[i*dof+j]-x[(i-1)*dof+j])/(x[(i-1)*dof+j]-x[(i-2)*dof+j]);
491: min[j] = PetscMin(r[j],2.0);
492: u[j] = x[(i-1)*dof+j]+PetscMax(0,PetscMin(min[j],alpha[0]+gamma[0]*r[j]))*(x[(i-1)*dof+j]-x[(i-2)*dof+j]);
493: }
494: (*ctx->physics.flux)(ctx->physics.user,u,ctx->flux,&maxspeed);
495: if (i > xs) {
496: for (j=0; j<dof; j++) f[(ifast-1)*dof+j] -= ctx->flux[j]/hf;
497: }
498: if (i < xs+xm) {
499: for (j=0; j<dof; j++) f[ifast*dof+j] += ctx->flux[j]/hf;
500: ifast++;
501: }
502: } else if (i == fs) {
503: u = &ctx->u[0];
504: alpha[0] = hf*hs/(hf+hf)/(hf+hf+hs);
505: gamma[0] = hf*(hf+hf)/(hf+hs)/(hf+hf+hs);
506: for (j=0; j<dof; j++) {
507: r[j] = (x[i*dof+j]-x[(i-1)*dof+j])/(x[(i-1)*dof+j]-x[(i-2)*dof+j]);
508: min[j] = PetscMin(r[j],2.0);
509: u[j] = x[(i-1)*dof+j]+PetscMax(0,PetscMin(min[j],alpha[0]+gamma[0]*r[j]))*(x[(i-1)*dof+j]-x[(i-2)*dof+j]);
510: }
511: (*ctx->physics.flux)(ctx->physics.user,u,ctx->flux,&maxspeed);
512: if (i > xs) {
513: for (j=0; j<dof; j++) f[(ifast-1)*dof+j] -= ctx->flux[j]/hf;
514: }
515: }
516: }
517: DMDAVecRestoreArray(da,Xloc,&x);
518: VecRestoreArray(F,&f);
519: DMRestoreLocalVector(da,&Xloc);
520: PetscFree4(r,min,alpha,gamma);
521: return(0);
522: }
524: /* --------------------------------- Finite Volume Solver for slow components ----------------------------------- */
526: PetscErrorCode FVSample(FVCtx *ctx,DM da,PetscReal time,Vec U)
527: {
529: PetscScalar *u,*uj,xj,xi;
530: PetscInt i,j,k,dof,xs,xm,Mx,count_slow,count_fast;
531: const PetscInt N=200;
534: if (!ctx->physics.sample) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Physics has not provided a sampling function");
535: DMDAGetInfo(da,0,&Mx,0,0,0,0,0,&dof,0,0,0,0,0);
536: DMDAGetCorners(da,&xs,0,0,&xm,0,0);
537: DMDAVecGetArray(da,U,&u);
538: PetscMalloc1(dof,&uj);
539: const PetscReal hs = (ctx->xmax-ctx->xmin)/2.0*(ctx->hratio+1.0)/Mx;
540: const PetscReal hf = (ctx->xmax-ctx->xmin)/2.0*(1.0+1.0/ctx->hratio)/Mx;
541: count_slow = Mx/(1+ctx->hratio);
542: count_fast = Mx-count_slow;
543: for (i=xs; i<xs+xm; i++) {
544: if (i*hs+0.5*hs<(ctx->xmax-ctx->xmin)*0.25) {
545: xi = ctx->xmin+0.5*hs+i*hs;
546: /* Integrate over cell i using trapezoid rule with N points. */
547: for (k=0; k<dof; k++) u[i*dof+k] = 0;
548: for (j=0; j<N+1; j++) {
549: xj = xi+hs*(j-N/2)/(PetscReal)N;
550: (*ctx->physics.sample)(ctx->physics.user,ctx->initial,ctx->bctype,ctx->xmin,ctx->xmax,time,xj,uj);
551: for (k=0; k<dof; k++) u[i*dof+k] += ((j==0 || j==N) ? 0.5 : 1.0)*uj[k]/N;
552: }
553: } else if ((ctx->xmax-ctx->xmin)*0.25+(i-count_slow/2)*hf+0.5*hf<(ctx->xmax-ctx->xmin)*0.75) {
554: xi = ctx->xmin+(ctx->xmax-ctx->xmin)*0.25+0.5*hf+(i-count_slow/2)*hf;
555: /* Integrate over cell i using trapezoid rule with N points. */
556: for (k=0; k<dof; k++) u[i*dof+k] = 0;
557: for (j=0; j<N+1; j++) {
558: xj = xi+hf*(j-N/2)/(PetscReal)N;
559: (*ctx->physics.sample)(ctx->physics.user,ctx->initial,ctx->bctype,ctx->xmin,ctx->xmax,time,xj,uj);
560: for (k=0; k<dof; k++) u[i*dof+k] += ((j==0 || j==N) ? 0.5 : 1.0)*uj[k]/N;
561: }
562: } else {
563: xi = ctx->xmin+(ctx->xmax-ctx->xmin)*0.75+0.5*hs+(i-count_slow/2-count_fast)*hs;
564: /* Integrate over cell i using trapezoid rule with N points. */
565: for (k=0; k<dof; k++) u[i*dof+k] = 0;
566: for (j=0; j<N+1; j++) {
567: xj = xi+hs*(j-N/2)/(PetscReal)N;
568: (*ctx->physics.sample)(ctx->physics.user,ctx->initial,ctx->bctype,ctx->xmin,ctx->xmax,time,xj,uj);
569: for (k=0; k<dof; k++) u[i*dof+k] += ((j==0 || j==N) ? 0.5 : 1.0)*uj[k]/N;
570: }
571: }
572: }
573: DMDAVecRestoreArray(da,U,&u);
574: PetscFree(uj);
575: return(0);
576: }
578: static PetscErrorCode SolutionStatsView(DM da,Vec X,PetscViewer viewer)
579: {
580: PetscErrorCode ierr;
581: PetscReal xmin,xmax;
582: PetscScalar sum,tvsum,tvgsum;
583: const PetscScalar *x;
584: PetscInt imin,imax,Mx,i,j,xs,xm,dof;
585: Vec Xloc;
586: PetscBool iascii;
589: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
590: if (iascii) {
591: /* PETSc lacks a function to compute total variation norm (difficult in multiple dimensions), we do it here */
592: DMGetLocalVector(da,&Xloc);
593: DMGlobalToLocalBegin(da,X,INSERT_VALUES,Xloc);
594: DMGlobalToLocalEnd(da,X,INSERT_VALUES,Xloc);
595: DMDAVecGetArrayRead(da,Xloc,(void*)&x);
596: DMDAGetCorners(da,&xs,0,0,&xm,0,0);
597: DMDAGetInfo(da,0,&Mx,0,0,0,0,0,&dof,0,0,0,0,0);
598: tvsum = 0;
599: for (i=xs; i<xs+xm; i++) {
600: for (j=0; j<dof; j++) tvsum += PetscAbsScalar(x[i*dof+j]-x[(i-1)*dof+j]);
601: }
602: MPI_Allreduce(&tvsum,&tvgsum,1,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)da));
603: DMDAVecRestoreArrayRead(da,Xloc,(void*)&x);
604: DMRestoreLocalVector(da,&Xloc);
606: VecMin(X,&imin,&xmin);
607: VecMax(X,&imax,&xmax);
608: VecSum(X,&sum);
609: PetscViewerASCIIPrintf(viewer,"Solution range [%g,%g] with minimum at %D, mean %g, ||x||_TV %g\n",(double)xmin,(double)xmax,imin,(double)(sum/Mx),(double)(tvgsum/Mx));
610: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Viewer type not supported");
611: return(0);
612: }
614: static PetscErrorCode SolutionErrorNorms(FVCtx *ctx,DM da,PetscReal t,Vec X,PetscReal *nrm1)
615: {
616: PetscErrorCode ierr;
617: Vec Y;
618: PetscInt i,Mx,count_slow=0,count_fast=0;
619: const PetscScalar *ptr_X,*ptr_Y;
622: VecGetSize(X,&Mx);
623: VecDuplicate(X,&Y);
624: FVSample(ctx,da,t,Y);
625: const PetscReal hs = (ctx->xmax-ctx->xmin)/2.0*(ctx->hratio+1.0)/Mx;
626: const PetscReal hf = (ctx->xmax-ctx->xmin)/2.0*(1.0+1.0/ctx->hratio)/Mx;
627: count_slow = (PetscReal)Mx/(1.0+ctx->hratio);
628: count_fast = Mx-count_slow;
629: VecGetArrayRead(X,&ptr_X);
630: VecGetArrayRead(Y,&ptr_Y);
631: for (i=0; i<Mx; i++) {
632: if (i < count_slow/2 || i > count_slow/2+count_fast-1) *nrm1 += hs*PetscAbs(ptr_X[i]-ptr_Y[i]);
633: else *nrm1 += hf*PetscAbs(ptr_X[i]-ptr_Y[i]);
634: }
635: VecRestoreArrayRead(X,&ptr_X);
636: VecRestoreArrayRead(Y,&ptr_Y);
637: VecDestroy(&Y);
638: return(0);
639: }
641: int main(int argc,char *argv[])
642: {
643: char physname[256] = "advect",final_fname[256] = "solution.m";
644: PetscFunctionList physics = 0;
645: MPI_Comm comm;
646: TS ts;
647: DM da;
648: Vec X,X0,R;
649: FVCtx ctx;
650: PetscInt i,k,dof,xs,xm,Mx,draw = 0,count_slow,count_fast,islow = 0,ifast = 0,*index_slow,*index_fast;
651: PetscBool view_final = PETSC_FALSE;
652: PetscReal ptime;
653: PetscErrorCode ierr;
655: PetscInitialize(&argc,&argv,0,help);if (ierr) return ierr;
656: comm = PETSC_COMM_WORLD;
657: PetscMemzero(&ctx,sizeof(ctx));
659: /* Register physical models to be available on the command line */
660: PetscFunctionListAdd(&physics,"advect",PhysicsCreate_Advect);
662: ctx.comm = comm;
663: ctx.cfl = 0.9;
664: ctx.bctype = FVBC_PERIODIC;
665: ctx.xmin = -1.0;
666: ctx.xmax = 1.0;
667: PetscOptionsBegin(comm,NULL,"Finite Volume solver options","");
668: PetscOptionsReal("-xmin","X min","",ctx.xmin,&ctx.xmin,NULL);
669: PetscOptionsReal("-xmax","X max","",ctx.xmax,&ctx.xmax,NULL);
670: PetscOptionsInt("-draw","Draw solution vector, bitwise OR of (1=initial,2=final,4=final error)","",draw,&draw,NULL);
671: PetscOptionsString("-view_final","Write final solution in ASCII MATLAB format to given file name","",final_fname,final_fname,sizeof(final_fname),&view_final);
672: PetscOptionsInt("-initial","Initial condition (depends on the physics)","",ctx.initial,&ctx.initial,NULL);
673: PetscOptionsBool("-exact","Compare errors with exact solution","",ctx.exact,&ctx.exact,NULL);
674: PetscOptionsBool("-simulation","Compare errors with reference solution","",ctx.simulation,&ctx.simulation,NULL);
675: PetscOptionsReal("-cfl","CFL number to time step at","",ctx.cfl,&ctx.cfl,NULL);
676: PetscOptionsEnum("-bc_type","Boundary condition","",FVBCTypes,(PetscEnum)ctx.bctype,(PetscEnum*)&ctx.bctype,NULL);
677: PetscOptionsInt("-hratio","Spacing ratio","",ctx.hratio,&ctx.hratio,NULL);
678: PetscOptionsEnd();
680: /* Choose the physics from the list of registered models */
681: {
682: PetscErrorCode (*r)(FVCtx*);
683: PetscFunctionListFind(physics,physname,&r);
684: if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Physics '%s' not found",physname);
685: /* Create the physics, will set the number of fields and their names */
686: (*r)(&ctx);
687: }
689: /* Create a DMDA to manage the parallel grid */
690: DMDACreate1d(comm,DM_BOUNDARY_PERIODIC,50,ctx.physics.dof,2,NULL,&da);
691: DMSetFromOptions(da);
692: DMSetUp(da);
693: /* Inform the DMDA of the field names provided by the physics. */
694: /* The names will be shown in the title bars when run with -ts_monitor_draw_solution */
695: for (i=0; i<ctx.physics.dof; i++) {
696: DMDASetFieldName(da,i,ctx.physics.fieldname[i]);
697: }
698: DMDAGetInfo(da,0,&Mx,0,0,0,0,0,&dof,0,0,0,0,0);
699: DMDAGetCorners(da,&xs,0,0,&xm,0,0);
701: /* Set coordinates of cell centers */
702: DMDASetUniformCoordinates(da,ctx.xmin+0.5*(ctx.xmax-ctx.xmin)/Mx,ctx.xmax+0.5*(ctx.xmax-ctx.xmin)/Mx,0,0,0,0);
704: /* Allocate work space for the Finite Volume solver (so it doesn't have to be reallocated on each function evaluation) */
705: PetscMalloc3(dof,&ctx.u,dof,&ctx.flux,dof,&ctx.speeds);
707: /* Create a vector to store the solution and to save the initial state */
708: DMCreateGlobalVector(da,&X);
709: VecDuplicate(X,&X0);
710: VecDuplicate(X,&R);
712: /* create index for slow parts and fast parts*/
713: count_slow = Mx/(1+ctx.hratio);
714: if (count_slow%2) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_USER,"Please adjust grid size Mx (-da_grid_x) and hratio (-hratio) so that Mx/(1+hartio) is even");
715: count_fast = Mx-count_slow;
716: ctx.sf = count_slow/2;
717: ctx.fs = ctx.sf + count_fast;
718: PetscMalloc1(xm*dof,&index_slow);
719: PetscMalloc1(xm*dof,&index_fast);
720: for (i=xs; i<xs+xm; i++) {
721: if (i < count_slow/2 || i > count_slow/2+count_fast-1)
722: for (k=0; k<dof; k++) index_slow[islow++] = i*dof+k;
723: else
724: for (k=0; k<dof; k++) index_fast[ifast++] = i*dof+k;
725: }
726: ISCreateGeneral(PETSC_COMM_WORLD,islow,index_slow,PETSC_COPY_VALUES,&ctx.iss);
727: ISCreateGeneral(PETSC_COMM_WORLD,ifast,index_fast,PETSC_COPY_VALUES,&ctx.isf);
729: /* Create a time-stepping object */
730: TSCreate(comm,&ts);
731: TSSetDM(ts,da);
732: TSSetRHSFunction(ts,R,FVRHSFunction,&ctx);
733: TSRHSSplitSetIS(ts,"slow",ctx.iss);
734: TSRHSSplitSetIS(ts,"fast",ctx.isf);
735: TSRHSSplitSetRHSFunction(ts,"slow",NULL,FVRHSFunctionslow,&ctx);
736: TSRHSSplitSetRHSFunction(ts,"fast",NULL,FVRHSFunctionfast,&ctx);
738: TSSetType(ts,TSMPRK);
739: TSSetMaxTime(ts,10);
740: TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);
742: /* Compute initial conditions and starting time step */
743: FVSample(&ctx,da,0,X0);
744: FVRHSFunction(ts,0,X0,X,(void*)&ctx); /* Initial function evaluation, only used to determine max speed */
745: VecCopy(X0,X); /* The function value was not used so we set X=X0 again */
746: TSSetTimeStep(ts,ctx.cfl/ctx.cfl_idt);
747: TSSetFromOptions(ts); /* Take runtime options */
748: SolutionStatsView(da,X,PETSC_VIEWER_STDOUT_WORLD);
749: {
750: PetscInt steps;
751: PetscScalar mass_initial,mass_final,mass_difference,mass_differenceg;
752: const PetscScalar *ptr_X,*ptr_X0;
753: const PetscReal hs = (ctx.xmax-ctx.xmin)/2.0/count_slow;
754: const PetscReal hf = (ctx.xmax-ctx.xmin)/2.0/count_fast;
755: TSSolve(ts,X);
756: TSGetSolveTime(ts,&ptime);
757: TSGetStepNumber(ts,&steps);
758: /* calculate the total mass at initial time and final time */
759: mass_initial = 0.0;
760: mass_final = 0.0;
761: DMDAVecGetArrayRead(da,X0,(void*)&ptr_X0);
762: DMDAVecGetArrayRead(da,X,(void*)&ptr_X);
763: for (i=xs; i<xs+xm; i++) {
764: if (i < ctx.sf || i > ctx.fs-1) {
765: for (k=0; k<dof; k++) {
766: mass_initial = mass_initial+hs*ptr_X0[i*dof+k];
767: mass_final = mass_final+hs*ptr_X[i*dof+k];
768: }
769: } else {
770: for (k=0; k<dof; k++) {
771: mass_initial = mass_initial+hf*ptr_X0[i*dof+k];
772: mass_final = mass_final+hf*ptr_X[i*dof+k];
773: }
774: }
775: }
776: DMDAVecRestoreArrayRead(da,X0,(void*)&ptr_X0);
777: DMDAVecRestoreArrayRead(da,X,(void*)&ptr_X);
778: mass_difference = mass_final-mass_initial;
779: MPI_Allreduce(&mass_difference,&mass_differenceg,1,MPIU_SCALAR,MPIU_SUM,comm);
780: PetscPrintf(comm,"Mass difference %g\n",(double)mass_differenceg);
781: PetscPrintf(comm,"Final time %g, steps %D\n",(double)ptime,steps);
782: if (ctx.exact) {
783: PetscReal nrm1 = 0;
784: SolutionErrorNorms(&ctx,da,ptime,X,&nrm1);
785: PetscPrintf(comm,"Error ||x-x_e||_1 %g\n",(double)nrm1);
786: }
787: if (ctx.simulation) {
788: PetscReal nrm1 = 0;
789: PetscViewer fd;
790: char filename[PETSC_MAX_PATH_LEN] = "binaryoutput";
791: Vec XR;
792: PetscBool flg;
793: const PetscScalar *ptr_XR;
794: PetscOptionsGetString(NULL,NULL,"-f",filename,sizeof(filename),&flg);
795: if (!flg) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_USER,"Must indicate binary file with the -f option");
796: PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_READ,&fd);
797: VecDuplicate(X0,&XR);
798: VecLoad(XR,fd);
799: PetscViewerDestroy(&fd);
800: VecGetArrayRead(X,&ptr_X);
801: VecGetArrayRead(XR,&ptr_XR);
802: for (i=0; i<Mx; i++) {
803: if (i < count_slow/2 || i > count_slow/2+count_fast-1) nrm1 = nrm1 + hs*PetscAbs(ptr_X[i]-ptr_XR[i]);
804: else nrm1 = nrm1 + hf*PetscAbs(ptr_X[i]-ptr_XR[i]);
805: }
806: VecRestoreArrayRead(X,&ptr_X);
807: VecRestoreArrayRead(XR,&ptr_XR);
808: PetscPrintf(comm,"Error ||x-x_e||_1 %g\n",(double)nrm1);
809: VecDestroy(&XR);
810: }
811: }
813: SolutionStatsView(da,X,PETSC_VIEWER_STDOUT_WORLD);
814: if (draw & 0x1) { VecView(X0,PETSC_VIEWER_DRAW_WORLD); }
815: if (draw & 0x2) { VecView(X,PETSC_VIEWER_DRAW_WORLD); }
816: if (draw & 0x4) {
817: Vec Y;
818: VecDuplicate(X,&Y);
819: FVSample(&ctx,da,ptime,Y);
820: VecAYPX(Y,-1,X);
821: VecView(Y,PETSC_VIEWER_DRAW_WORLD);
822: VecDestroy(&Y);
823: }
825: if (view_final) {
826: PetscViewer viewer;
827: PetscViewerASCIIOpen(PETSC_COMM_WORLD,final_fname,&viewer);
828: PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_MATLAB);
829: VecView(X,viewer);
830: PetscViewerPopFormat(viewer);
831: PetscViewerDestroy(&viewer);
832: }
834: /* Clean up */
835: (*ctx.physics.destroy)(ctx.physics.user);
836: for (i=0; i<ctx.physics.dof; i++) {PetscFree(ctx.physics.fieldname[i]);}
837: PetscFree3(ctx.u,ctx.flux,ctx.speeds);
838: ISDestroy(&ctx.iss);
839: ISDestroy(&ctx.isf);
840: VecDestroy(&X);
841: VecDestroy(&X0);
842: VecDestroy(&R);
843: DMDestroy(&da);
844: TSDestroy(&ts);
845: PetscFree(index_slow);
846: PetscFree(index_fast);
847: PetscFunctionListDestroy(&physics);
848: PetscFinalize();
849: return ierr;
850: }
852: /*TEST
854: build:
855: requires: !complex
857: test:
858: args: -da_grid_x 60 -initial 7 -xmin -1 -xmax 1 -hratio 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
860: test:
861: suffix: 2
862: args: -da_grid_x 60 -initial 7 -xmin -1 -xmax 1 -hratio 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
863: output_file: output/ex7_1.out
865: test:
866: suffix: 3
867: args: -da_grid_x 60 -initial 7 -xmin -1 -xmax 1 -hratio 2 -ts_dt 0.025 -ts_max_steps 24 -ts_type mprk -ts_mprk_type 2a22 -ts_use_splitrhsfunction 0
869: test:
870: suffix: 4
871: args: -da_grid_x 60 -initial 7 -xmin -1 -xmax 1 -hratio 2 -ts_dt 0.025 -ts_max_steps 24 -ts_type mprk -ts_mprk_type 2a22 -ts_use_splitrhsfunction 1
872: output_file: output/ex7_3.out
874: test:
875: suffix: 5
876: nsize: 2
877: args: -da_grid_x 60 -initial 7 -xmin -1 -xmax 1 -hratio 2 -ts_dt 0.025 -ts_max_steps 24 -ts_type mprk -ts_mprk_type 2a22 -ts_use_splitrhsfunction 1
878: output_file: output/ex7_3.out
879: TEST*/