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 coarse\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: 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: {
71: PetscFree(vctx);
72: return 0;
73: }
75: /* --------------------------------- Advection ----------------------------------- */
76: typedef struct {
77: PetscReal a; /* advective velocity */
78: } AdvectCtx;
80: static PetscErrorCode PhysicsFlux_Advect(void *vctx,const PetscScalar *u,PetscScalar *flux,PetscReal *maxspeed)
81: {
82: AdvectCtx *ctx = (AdvectCtx*)vctx;
83: PetscReal speed;
86: speed = ctx->a;
87: flux[0] = speed*u[0];
88: *maxspeed = speed;
89: return 0;
90: }
92: static PetscErrorCode PhysicsSample_Advect(void *vctx,PetscInt initial,FVBCType bctype,PetscReal xmin,PetscReal xmax,PetscReal t,PetscReal x,PetscReal *u)
93: {
94: AdvectCtx *ctx = (AdvectCtx*)vctx;
95: PetscReal a = ctx->a,x0;
98: switch (bctype) {
99: case FVBC_OUTFLOW: x0 = x-a*t; break;
100: case FVBC_PERIODIC: x0 = RangeMod(x-a*t,xmin,xmax); break;
101: default: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"unknown BCType");
102: }
103: switch (initial) {
104: case 0: u[0] = (x0 < 0) ? 1 : -1; break;
105: case 1: u[0] = (x0 < 0) ? -1 : 1; break;
106: case 2: u[0] = (0 < x0 && x0 < 1) ? 1 : 0; break;
107: case 3: u[0] = PetscSinReal(2*PETSC_PI*x0); break;
108: case 4: u[0] = PetscAbs(x0); break;
109: case 5: u[0] = (x0 < 0 || x0 > 0.5) ? 0 : PetscSqr(PetscSinReal(2*PETSC_PI*x0)); break;
110: case 6: u[0] = (x0 < 0) ? 0 : ((x0 < 1) ? x0 : ((x0 < 2) ? 2-x0 : 0)); break;
111: case 7: u[0] = PetscPowReal(PetscSinReal(PETSC_PI*x0),10.0);break;
112: default: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"unknown initial condition");
113: }
114: return 0;
115: }
117: static PetscErrorCode PhysicsCreate_Advect(FVCtx *ctx)
118: {
120: AdvectCtx *user;
123: PetscNew(&user);
124: ctx->physics.sample = PhysicsSample_Advect;
125: ctx->physics.flux = PhysicsFlux_Advect;
126: ctx->physics.destroy = PhysicsDestroy_SimpleFree;
127: ctx->physics.user = user;
128: ctx->physics.dof = 1;
129: PetscStrallocpy("u",&ctx->physics.fieldname[0]);
130: user->a = 1;
131: PetscOptionsBegin(ctx->comm,ctx->prefix,"Options for advection","");
132: {
133: PetscOptionsReal("-physics_advect_a","Speed","",user->a,&user->a,NULL);
134: }
135: PetscOptionsEnd();
136: return 0;
137: }
139: /* --------------------------------- Finite Volume Solver ----------------------------------- */
141: static PetscErrorCode FVRHSFunction(TS ts,PetscReal time,Vec X,Vec F,void *vctx)
142: {
143: FVCtx *ctx = (FVCtx*)vctx;
144: PetscInt i,j,Mx,dof,xs,xm,sf = ctx->sf,fs = ctx->fs;
145: PetscReal hf,hs,cfl_idt = 0;
146: PetscScalar *x,*f,*r,*min,*alpha,*gamma;
147: Vec Xloc;
148: DM da;
151: TSGetDM(ts,&da);
152: DMGetLocalVector(da,&Xloc); /* Xloc contains ghost points */
153: DMDAGetInfo(da,0,&Mx,0,0,0,0,0,&dof,0,0,0,0,0); /* Mx is the number of center points */
154: hs = (ctx->xmax-ctx->xmin)/2.0*(ctx->hratio+1.0)/Mx;
155: hf = (ctx->xmax-ctx->xmin)/2.0*(1.0+1.0/ctx->hratio)/Mx;
156: DMGlobalToLocalBegin(da,X,INSERT_VALUES,Xloc); /* X is solution vector which does not contain ghost points */
157: DMGlobalToLocalEnd(da,X,INSERT_VALUES,Xloc);
158: VecZeroEntries(F); /* F is the right hand side function corresponds to center points */
159: DMDAVecGetArray(da,Xloc,&x);
160: DMDAVecGetArray(da,F,&f);
161: DMDAGetCorners(da,&xs,0,0,&xm,0,0);
162: PetscMalloc4(dof,&r,dof,&min,dof,&alpha,dof,&gamma);
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: }
173: for (i=xs; i<xs+xm+1; i++) {
174: PetscReal maxspeed;
175: PetscScalar *u;
176: if (i < sf || i > fs+1) {
177: u = &ctx->u[0];
178: alpha[0] = 1.0/6.0;
179: gamma[0] = 1.0/3.0;
180: for (j=0; j<dof; j++) {
181: r[j] = (x[i*dof+j]-x[(i-1)*dof+j])/(x[(i-1)*dof+j]-x[(i-2)*dof+j]);
182: min[j] = PetscMin(r[j],2.0);
183: 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]);
184: }
185: (*ctx->physics.flux)(ctx->physics.user,u,ctx->flux,&maxspeed);
186: cfl_idt = PetscMax(cfl_idt,PetscAbsScalar(maxspeed/hs));
187: if (i > xs) {
188: for (j=0; j<dof; j++) f[(i-1)*dof+j] -= ctx->flux[j]/hs;
189: }
190: if (i < xs+xm) {
191: for (j=0; j<dof; j++) f[i*dof+j] += ctx->flux[j]/hs;
192: }
193: } else if (i == sf) {
194: u = &ctx->u[0];
195: alpha[0] = hs*hf/(hs+hs)/(hs+hs+hf);
196: gamma[0] = hs*(hs+hs)/(hs+hf)/(hs+hs+hf);
197: for (j=0; j<dof; j++) {
198: r[j] = (x[i*dof+j]-x[(i-1)*dof+j])/(x[(i-1)*dof+j]-x[(i-2)*dof+j]);
199: min[j] = PetscMin(r[j],2.0);
200: 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]);
201: }
202: (*ctx->physics.flux)(ctx->physics.user,u,ctx->flux,&maxspeed);
203: if (i > xs) {
204: for (j=0; j<dof; j++) f[(i-1)*dof+j] -= ctx->flux[j]/hs;
205: }
206: if (i < xs+xm) {
207: for (j=0; j<dof; j++) f[i*dof+j] += ctx->flux[j]/hf;
208: }
209: } else if (i == sf+1) {
210: u = &ctx->u[0];
211: alpha[0] = hf*hf/(hs+hf)/(hs+hf+hf);
212: gamma[0] = hf*(hs+hf)/(hf+hf)/(hs+hf+hf);
213: for (j=0; j<dof; j++) {
214: r[j] = (x[i*dof+j] - x[(i-1)*dof+j])/(x[(i-1)*dof+j]-x[(i-2)*dof+j]);
215: min[j] = PetscMin(r[j],2.0);
216: 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]);
217: }
218: (*ctx->physics.flux)(ctx->physics.user,u,ctx->flux,&maxspeed);
219: if (i > xs) {
220: for (j=0; j<dof; j++) f[(i-1)*dof+j] -= ctx->flux[j]/hf;
221: }
222: if (i < xs+xm) {
223: for (j=0; j<dof; j++) f[i*dof+j] += ctx->flux[j]/hf;
224: }
225: } else if (i > sf+1 && i < fs) {
226: u = &ctx->u[0];
227: alpha[0] = 1.0/6.0;
228: gamma[0] = 1.0/3.0;
229: for (j=0; j<dof; j++) {
230: r[j] = (x[i*dof+j] - x[(i-1)*dof+j])/(x[(i-1)*dof+j]-x[(i-2)*dof+j]);
231: min[j] = PetscMin(r[j],2.0);
232: 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]);
233: }
234: (*ctx->physics.flux)(ctx->physics.user,u,ctx->flux,&maxspeed);
235: if (i > xs) {
236: for (j=0; j<dof; j++) f[(i-1)*dof+j] -= ctx->flux[j]/hf;
237: }
238: if (i < xs+xm) {
239: for (j=0; j<dof; j++) f[i*dof+j] += ctx->flux[j]/hf;
240: }
241: } else if (i == fs) {
242: u = &ctx->u[0];
243: alpha[0] = hf*hs/(hf+hf)/(hf+hf+hs);
244: gamma[0] = hf*(hf+hf)/(hf+hs)/(hf+hf+hs);
245: for (j=0; j<dof; j++) {
246: r[j] = (x[i*dof+j] - x[(i-1)*dof+j])/(x[(i-1)*dof+j]-x[(i-2)*dof+j]);
247: min[j] = PetscMin(r[j],2.0);
248: 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]);
249: }
250: (*ctx->physics.flux)(ctx->physics.user,u,ctx->flux,&maxspeed);
251: if (i > xs) {
252: for (j=0; j<dof; j++) f[(i-1)*dof+j] -= ctx->flux[j]/hf;
253: }
254: if (i < xs+xm) {
255: for (j=0; j<dof; j++) f[i*dof+j] += ctx->flux[j]/hs;
256: }
257: } else if (i == fs+1) {
258: u = &ctx->u[0];
259: alpha[0] = hs*hs/(hf+hs)/(hf+hs+hs);
260: gamma[0] = hs*(hf+hs)/(hs+hs)/(hf+hs+hs);
261: for (j=0; j<dof; j++) {
262: r[j] = (x[i*dof+j] - x[(i-1)*dof+j])/(x[(i-1)*dof+j]-x[(i-2)*dof+j]);
263: min[j] = PetscMin(r[j],2.0);
264: 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]);
265: }
266: (*ctx->physics.flux)(ctx->physics.user,u,ctx->flux,&maxspeed);
267: if (i > xs) {
268: for (j=0; j<dof; j++) f[(i-1)*dof+j] -= ctx->flux[j]/hs;
269: }
270: if (i < xs+xm) {
271: for (j=0; j<dof; j++) f[i*dof+j] += ctx->flux[j]/hs;
272: }
273: }
274: }
275: DMDAVecRestoreArray(da,Xloc,&x);
276: DMDAVecRestoreArray(da,F,&f);
277: DMRestoreLocalVector(da,&Xloc);
278: MPI_Allreduce(&cfl_idt,&ctx->cfl_idt,1,MPIU_SCALAR,MPIU_MAX,PetscObjectComm((PetscObject)da));
279: if (0) {
280: /* We need a way to inform the TS of a CFL constraint, this is a debugging fragment */
281: PetscReal dt,tnow;
282: TSGetTimeStep(ts,&dt);
283: TSGetTime(ts,&tnow);
284: if (dt > 0.5/ctx->cfl_idt) {
285: PetscPrintf(ctx->comm,"Stability constraint exceeded at t=%g, dt %g > %g\n",(double)tnow,(double)dt,(double)(0.5/ctx->cfl_idt));
286: }
287: }
288: PetscFree4(r,min,alpha,gamma);
289: return 0;
290: }
292: static PetscErrorCode FVRHSFunctionslow(TS ts,PetscReal time,Vec X,Vec F,void *vctx)
293: {
294: FVCtx *ctx = (FVCtx*)vctx;
295: PetscInt i,j,Mx,dof,xs,xm,islow = 0,sf = ctx->sf,fs = ctx->fs;
296: PetscReal hf,hs;
297: PetscScalar *x,*f,*r,*min,*alpha,*gamma;
298: Vec Xloc;
299: DM da;
302: TSGetDM(ts,&da);
303: DMGetLocalVector(da,&Xloc); /* Xloc contains ghost points */
304: DMDAGetInfo(da,0,&Mx,0,0,0,0,0,&dof,0,0,0,0,0); /* Mx is the number of center points */
305: hs = (ctx->xmax-ctx->xmin)/2.0*(ctx->hratio+1.0)/Mx;
306: hf = (ctx->xmax-ctx->xmin)/2.0*(1.0+1.0/ctx->hratio)/Mx;
307: DMGlobalToLocalBegin(da,X,INSERT_VALUES,Xloc); /* X is solution vector which does not contain ghost points */
308: DMGlobalToLocalEnd (da,X,INSERT_VALUES,Xloc);
309: VecZeroEntries(F); /* F is the right hand side function corresponds to center points */
310: DMDAVecGetArray(da,Xloc,&x);
311: VecGetArray(F,&f);
312: DMDAGetCorners(da,&xs,0,0,&xm,0,0);
313: PetscMalloc4(dof,&r,dof,&min,dof,&alpha,dof,&gamma);
315: if (ctx->bctype == FVBC_OUTFLOW) {
316: for (i=xs-2; i<0; i++) {
317: for (j=0; j<dof; j++) x[i*dof+j] = x[j];
318: }
319: for (i=Mx; i<xs+xm+2; i++) {
320: for (j=0; j<dof; j++) x[i*dof+j] = x[(xs+xm-1)*dof+j];
321: }
322: }
324: for (i=xs; i<xs+xm+1; i++) {
325: PetscReal maxspeed;
326: PetscScalar *u;
327: if (i < sf) {
328: u = &ctx->u[0];
329: alpha[0] = 1.0/6.0;
330: gamma[0] = 1.0/3.0;
331: for (j=0; j<dof; j++) {
332: r[j] = (x[i*dof+j] - x[(i-1)*dof+j])/(x[(i-1)*dof+j]-x[(i-2)*dof+j]);
333: min[j] = PetscMin(r[j],2.0);
334: 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]);
335: }
336: (*ctx->physics.flux)(ctx->physics.user,u,ctx->flux,&maxspeed);
337: if (i > xs) {
338: for (j=0; j<dof; j++) f[(islow-1)*dof+j] -= ctx->flux[j]/hs;
339: }
340: if (i < xs+xm) {
341: for (j=0; j<dof; j++) f[islow*dof+j] += ctx->flux[j]/hs;
342: islow++;
343: }
344: } else if (i == sf) {
345: u = &ctx->u[0];
346: alpha[0] = hs*hf/(hs+hs)/(hs+hs+hf);
347: gamma[0] = hs*(hs+hs)/(hs+hf)/(hs+hs+hf);
348: for (j=0; j<dof; j++) {
349: r[j] = (x[i*dof+j] - x[(i-1)*dof+j])/(x[(i-1)*dof+j]-x[(i-2)*dof+j]);
350: min[j] = PetscMin(r[j],2.0);
351: 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]);
352: }
353: (*ctx->physics.flux)(ctx->physics.user,u,ctx->flux,&maxspeed);
354: if (i > xs) {
355: for (j=0; j<dof; j++) f[(islow-1)*dof+j] -= ctx->flux[j]/hs;
356: }
357: } else if (i == fs) {
358: u = &ctx->u[0];
359: alpha[0] = hf*hs/(hf+hf)/(hf+hf+hs);
360: gamma[0] = hf*(hf+hf)/(hf+hs)/(hf+hf+hs);
361: for (j=0; j<dof; j++) {
362: r[j] = (x[i*dof+j] - x[(i-1)*dof+j])/(x[(i-1)*dof+j]-x[(i-2)*dof+j]);
363: min[j] = PetscMin(r[j],2.0);
364: 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]);
365: }
366: (*ctx->physics.flux)(ctx->physics.user,u,ctx->flux,&maxspeed);
367: if (i < xs+xm) {
368: for (j=0; j<dof; j++) f[islow*dof+j] += ctx->flux[j]/hs;
369: islow++;
370: }
371: } else if (i == fs+1) {
372: u = &ctx->u[0];
373: alpha[0] = hs*hs/(hf+hs)/(hf+hs+hs);
374: gamma[0] = hs*(hf+hs)/(hs+hs)/(hf+hs+hs);
375: for (j=0; j<dof; j++) {
376: r[j] = (x[i*dof+j] - x[(i-1)*dof+j])/(x[(i-1)*dof+j]-x[(i-2)*dof+j]);
377: min[j] = PetscMin(r[j],2.0);
378: 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]);
379: }
380: (*ctx->physics.flux)(ctx->physics.user,u,ctx->flux,&maxspeed);
381: if (i > xs) {
382: for (j=0; j<dof; j++) f[(islow-1)*dof+j] -= ctx->flux[j]/hs;
383: }
384: if (i < xs+xm) {
385: for (j=0; j<dof; j++) f[islow*dof+j] += ctx->flux[j]/hs;
386: islow++;
387: }
388: } else if (i > fs+1) {
389: u = &ctx->u[0];
390: alpha[0] = 1.0/6.0;
391: gamma[0] = 1.0/3.0;
392: for (j=0; j<dof; j++) {
393: r[j] = (x[i*dof+j] - x[(i-1)*dof+j])/(x[(i-1)*dof+j]-x[(i-2)*dof+j]);
394: min[j] = PetscMin(r[j],2.0);
395: 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]);
396: }
397: (*ctx->physics.flux)(ctx->physics.user,u,ctx->flux,&maxspeed);
398: if (i > xs) {
399: for (j=0; j<dof; j++) f[(islow-1)*dof+j] -= ctx->flux[j]/hs;
400: }
401: if (i < xs+xm) {
402: for (j=0; j<dof; j++) f[islow*dof+j] += ctx->flux[j]/hs;
403: islow++;
404: }
405: }
406: }
407: DMDAVecRestoreArray(da,Xloc,&x);
408: VecRestoreArray(F,&f);
409: DMRestoreLocalVector(da,&Xloc);
410: PetscFree4(r,min,alpha,gamma);
411: return 0;
412: }
414: static PetscErrorCode FVRHSFunctionfast(TS ts,PetscReal time,Vec X,Vec F,void *vctx)
415: {
416: FVCtx *ctx = (FVCtx*)vctx;
417: PetscInt i,j,Mx,dof,xs,xm,ifast = 0,sf = ctx->sf,fs = ctx->fs;
418: PetscReal hf,hs;
419: PetscScalar *x,*f,*r,*min,*alpha,*gamma;
420: Vec Xloc;
421: DM da;
424: TSGetDM(ts,&da);
425: DMGetLocalVector(da,&Xloc); /* Xloc contains ghost points */
426: DMDAGetInfo(da,0,&Mx,0,0,0,0,0,&dof,0,0,0,0,0); /* Mx is the number of center points */
427: hs = (ctx->xmax-ctx->xmin)/2.0*(ctx->hratio+1.0)/Mx;
428: hf = (ctx->xmax-ctx->xmin)/2.0*(1.0+1.0/ctx->hratio)/Mx;
429: DMGlobalToLocalBegin(da,X,INSERT_VALUES,Xloc); /* X is solution vector which does not contain ghost points */
430: DMGlobalToLocalEnd(da,X,INSERT_VALUES,Xloc);
431: VecZeroEntries(F); /* F is the right hand side function corresponds to center points */
432: DMDAVecGetArray(da,Xloc,&x);
433: VecGetArray(F,&f);
434: DMDAGetCorners(da,&xs,0,0,&xm,0,0);
435: PetscMalloc4(dof,&r,dof,&min,dof,&alpha,dof,&gamma);
437: if (ctx->bctype == FVBC_OUTFLOW) {
438: for (i=xs-2; i<0; i++) {
439: for (j=0; j<dof; j++) x[i*dof+j] = x[j];
440: }
441: for (i=Mx; i<xs+xm+2; i++) {
442: for (j=0; j<dof; j++) x[i*dof+j] = x[(xs+xm-1)*dof+j];
443: }
444: }
446: for (i=xs; i<xs+xm+1; i++) {
447: PetscReal maxspeed;
448: PetscScalar *u;
449: if (i == sf) {
450: u = &ctx->u[0];
451: alpha[0] = hs*hf/(hs+hs)/(hs+hs+hf);
452: gamma[0] = hs*(hs+hs)/(hs+hf)/(hs+hs+hf);
453: for (j=0; j<dof; j++) {
454: r[j] = (x[i*dof+j]-x[(i-1)*dof+j])/(x[(i-1)*dof+j]-x[(i-2)*dof+j]);
455: min[j] = PetscMin(r[j],2.0);
456: 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]);
457: }
458: (*ctx->physics.flux)(ctx->physics.user,u,ctx->flux,&maxspeed);
459: if (i < xs+xm) {
460: for (j=0; j<dof; j++) f[ifast*dof+j] += ctx->flux[j]/hf;
461: ifast++;
462: }
463: } else if (i == sf+1) {
464: u = &ctx->u[0];
465: alpha[0] = hf*hf/(hs+hf)/(hs+hf+hf);
466: gamma[0] = hf*(hs+hf)/(hf+hf)/(hs+hf+hf);
467: for (j=0; j<dof; j++) {
468: r[j] = (x[i*dof+j]-x[(i-1)*dof+j])/(x[(i-1)*dof+j]-x[(i-2)*dof+j]);
469: min[j] = PetscMin(r[j],2.0);
470: 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]);
471: }
472: (*ctx->physics.flux)(ctx->physics.user,u,ctx->flux,&maxspeed);
473: if (i > xs) {
474: for (j=0; j<dof; j++) f[(ifast-1)*dof+j] -= ctx->flux[j]/hf;
475: }
476: if (i < xs+xm) {
477: for (j=0; j<dof; j++) f[ifast*dof+j] += ctx->flux[j]/hf;
478: ifast++;
479: }
480: } else if (i > sf+1 && i < fs) {
481: u = &ctx->u[0];
482: alpha[0] = 1.0/6.0;
483: gamma[0] = 1.0/3.0;
484: for (j=0; j<dof; j++) {
485: r[j] = (x[i*dof+j]-x[(i-1)*dof+j])/(x[(i-1)*dof+j]-x[(i-2)*dof+j]);
486: min[j] = PetscMin(r[j],2.0);
487: 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]);
488: }
489: (*ctx->physics.flux)(ctx->physics.user,u,ctx->flux,&maxspeed);
490: if (i > xs) {
491: for (j=0; j<dof; j++) f[(ifast-1)*dof+j] -= ctx->flux[j]/hf;
492: }
493: if (i < xs+xm) {
494: for (j=0; j<dof; j++) f[ifast*dof+j] += ctx->flux[j]/hf;
495: ifast++;
496: }
497: } else if (i == fs) {
498: u = &ctx->u[0];
499: alpha[0] = hf*hs/(hf+hf)/(hf+hf+hs);
500: gamma[0] = hf*(hf+hf)/(hf+hs)/(hf+hf+hs);
501: for (j=0; j<dof; j++) {
502: r[j] = (x[i*dof+j]-x[(i-1)*dof+j])/(x[(i-1)*dof+j]-x[(i-2)*dof+j]);
503: min[j] = PetscMin(r[j],2.0);
504: 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]);
505: }
506: (*ctx->physics.flux)(ctx->physics.user,u,ctx->flux,&maxspeed);
507: if (i > xs) {
508: for (j=0; j<dof; j++) f[(ifast-1)*dof+j] -= ctx->flux[j]/hf;
509: }
510: }
511: }
512: DMDAVecRestoreArray(da,Xloc,&x);
513: VecRestoreArray(F,&f);
514: DMRestoreLocalVector(da,&Xloc);
515: PetscFree4(r,min,alpha,gamma);
516: return 0;
517: }
519: /* --------------------------------- Finite Volume Solver for slow components ----------------------------------- */
521: PetscErrorCode FVSample(FVCtx *ctx,DM da,PetscReal time,Vec U)
522: {
523: PetscScalar *u,*uj,xj,xi;
524: PetscInt i,j,k,dof,xs,xm,Mx,count_slow,count_fast;
525: const PetscInt N=200;
529: DMDAGetInfo(da,0,&Mx,0,0,0,0,0,&dof,0,0,0,0,0);
530: DMDAGetCorners(da,&xs,0,0,&xm,0,0);
531: DMDAVecGetArray(da,U,&u);
532: PetscMalloc1(dof,&uj);
533: const PetscReal hs = (ctx->xmax-ctx->xmin)/2.0*(ctx->hratio+1.0)/Mx;
534: const PetscReal hf = (ctx->xmax-ctx->xmin)/2.0*(1.0+1.0/ctx->hratio)/Mx;
535: count_slow = Mx/(1+ctx->hratio);
536: count_fast = Mx-count_slow;
537: for (i=xs; i<xs+xm; i++) {
538: if (i*hs+0.5*hs<(ctx->xmax-ctx->xmin)*0.25) {
539: xi = ctx->xmin+0.5*hs+i*hs;
540: /* Integrate over cell i using trapezoid rule with N points. */
541: for (k=0; k<dof; k++) u[i*dof+k] = 0;
542: for (j=0; j<N+1; j++) {
543: xj = xi+hs*(j-N/2)/(PetscReal)N;
544: (*ctx->physics.sample)(ctx->physics.user,ctx->initial,ctx->bctype,ctx->xmin,ctx->xmax,time,xj,uj);
545: for (k=0; k<dof; k++) u[i*dof+k] += ((j==0 || j==N) ? 0.5 : 1.0)*uj[k]/N;
546: }
547: } else if ((ctx->xmax-ctx->xmin)*0.25+(i-count_slow/2)*hf+0.5*hf<(ctx->xmax-ctx->xmin)*0.75) {
548: xi = ctx->xmin+(ctx->xmax-ctx->xmin)*0.25+0.5*hf+(i-count_slow/2)*hf;
549: /* Integrate over cell i using trapezoid rule with N points. */
550: for (k=0; k<dof; k++) u[i*dof+k] = 0;
551: for (j=0; j<N+1; j++) {
552: xj = xi+hf*(j-N/2)/(PetscReal)N;
553: (*ctx->physics.sample)(ctx->physics.user,ctx->initial,ctx->bctype,ctx->xmin,ctx->xmax,time,xj,uj);
554: for (k=0; k<dof; k++) u[i*dof+k] += ((j==0 || j==N) ? 0.5 : 1.0)*uj[k]/N;
555: }
556: } else {
557: xi = ctx->xmin+(ctx->xmax-ctx->xmin)*0.75+0.5*hs+(i-count_slow/2-count_fast)*hs;
558: /* Integrate over cell i using trapezoid rule with N points. */
559: for (k=0; k<dof; k++) u[i*dof+k] = 0;
560: for (j=0; j<N+1; j++) {
561: xj = xi+hs*(j-N/2)/(PetscReal)N;
562: (*ctx->physics.sample)(ctx->physics.user,ctx->initial,ctx->bctype,ctx->xmin,ctx->xmax,time,xj,uj);
563: for (k=0; k<dof; k++) u[i*dof+k] += ((j==0 || j==N) ? 0.5 : 1.0)*uj[k]/N;
564: }
565: }
566: }
567: DMDAVecRestoreArray(da,U,&u);
568: PetscFree(uj);
569: return 0;
570: }
572: static PetscErrorCode SolutionStatsView(DM da,Vec X,PetscViewer viewer)
573: {
574: PetscReal xmin,xmax;
575: PetscScalar sum,tvsum,tvgsum;
576: const PetscScalar *x;
577: PetscInt imin,imax,Mx,i,j,xs,xm,dof;
578: Vec Xloc;
579: PetscBool iascii;
582: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
583: if (iascii) {
584: /* PETSc lacks a function to compute total variation norm (difficult in multiple dimensions), we do it here */
585: DMGetLocalVector(da,&Xloc);
586: DMGlobalToLocalBegin(da,X,INSERT_VALUES,Xloc);
587: DMGlobalToLocalEnd(da,X,INSERT_VALUES,Xloc);
588: DMDAVecGetArrayRead(da,Xloc,(void*)&x);
589: DMDAGetCorners(da,&xs,0,0,&xm,0,0);
590: DMDAGetInfo(da,0,&Mx,0,0,0,0,0,&dof,0,0,0,0,0);
591: tvsum = 0;
592: for (i=xs; i<xs+xm; i++) {
593: for (j=0; j<dof; j++) tvsum += PetscAbsScalar(x[i*dof+j]-x[(i-1)*dof+j]);
594: }
595: MPI_Allreduce(&tvsum,&tvgsum,1,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)da));
596: DMDAVecRestoreArrayRead(da,Xloc,(void*)&x);
597: DMRestoreLocalVector(da,&Xloc);
599: VecMin(X,&imin,&xmin);
600: VecMax(X,&imax,&xmax);
601: VecSum(X,&sum);
602: 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));
603: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Viewer type not supported");
604: return 0;
605: }
607: static PetscErrorCode SolutionErrorNorms(FVCtx *ctx,DM da,PetscReal t,Vec X,PetscReal *nrm1)
608: {
609: Vec Y;
610: PetscInt i,Mx,count_slow=0,count_fast=0;
611: const PetscScalar *ptr_X,*ptr_Y;
614: VecGetSize(X,&Mx);
615: VecDuplicate(X,&Y);
616: FVSample(ctx,da,t,Y);
617: const PetscReal hs = (ctx->xmax-ctx->xmin)/2.0*(ctx->hratio+1.0)/Mx;
618: const PetscReal hf = (ctx->xmax-ctx->xmin)/2.0*(1.0+1.0/ctx->hratio)/Mx;
619: count_slow = (PetscReal)Mx/(1.0+ctx->hratio);
620: count_fast = Mx-count_slow;
621: VecGetArrayRead(X,&ptr_X);
622: VecGetArrayRead(Y,&ptr_Y);
623: for (i=0; i<Mx; i++) {
624: if (i < count_slow/2 || i > count_slow/2+count_fast-1) *nrm1 += hs*PetscAbs(ptr_X[i]-ptr_Y[i]);
625: else *nrm1 += hf*PetscAbs(ptr_X[i]-ptr_Y[i]);
626: }
627: VecRestoreArrayRead(X,&ptr_X);
628: VecRestoreArrayRead(Y,&ptr_Y);
629: VecDestroy(&Y);
630: return 0;
631: }
633: int main(int argc,char *argv[])
634: {
635: char physname[256] = "advect",final_fname[256] = "solution.m";
636: PetscFunctionList physics = 0;
637: MPI_Comm comm;
638: TS ts;
639: DM da;
640: Vec X,X0,R;
641: FVCtx ctx;
642: PetscInt i,k,dof,xs,xm,Mx,draw = 0,count_slow,count_fast,islow = 0,ifast = 0,*index_slow,*index_fast;
643: PetscBool view_final = PETSC_FALSE;
644: PetscReal ptime;
645: PetscErrorCode ierr;
647: PetscInitialize(&argc,&argv,0,help);
648: comm = PETSC_COMM_WORLD;
649: PetscMemzero(&ctx,sizeof(ctx));
651: /* Register physical models to be available on the command line */
652: PetscFunctionListAdd(&physics,"advect",PhysicsCreate_Advect);
654: ctx.comm = comm;
655: ctx.cfl = 0.9;
656: ctx.bctype = FVBC_PERIODIC;
657: ctx.xmin = -1.0;
658: ctx.xmax = 1.0;
659: PetscOptionsBegin(comm,NULL,"Finite Volume solver options","");
660: PetscOptionsReal("-xmin","X min","",ctx.xmin,&ctx.xmin,NULL);
661: PetscOptionsReal("-xmax","X max","",ctx.xmax,&ctx.xmax,NULL);
662: PetscOptionsInt("-draw","Draw solution vector, bitwise OR of (1=initial,2=final,4=final error)","",draw,&draw,NULL);
663: PetscOptionsString("-view_final","Write final solution in ASCII MATLAB format to given file name","",final_fname,final_fname,sizeof(final_fname),&view_final);
664: PetscOptionsInt("-initial","Initial condition (depends on the physics)","",ctx.initial,&ctx.initial,NULL);
665: PetscOptionsBool("-exact","Compare errors with exact solution","",ctx.exact,&ctx.exact,NULL);
666: PetscOptionsBool("-simulation","Compare errors with reference solution","",ctx.simulation,&ctx.simulation,NULL);
667: PetscOptionsReal("-cfl","CFL number to time step at","",ctx.cfl,&ctx.cfl,NULL);
668: PetscOptionsEnum("-bc_type","Boundary condition","",FVBCTypes,(PetscEnum)ctx.bctype,(PetscEnum*)&ctx.bctype,NULL);
669: PetscOptionsInt("-hratio","Spacing ratio","",ctx.hratio,&ctx.hratio,NULL);
670: PetscOptionsEnd();
672: /* Choose the physics from the list of registered models */
673: {
674: PetscErrorCode (*r)(FVCtx*);
675: PetscFunctionListFind(physics,physname,&r);
677: /* Create the physics, will set the number of fields and their names */
678: (*r)(&ctx);
679: }
681: /* Create a DMDA to manage the parallel grid */
682: DMDACreate1d(comm,DM_BOUNDARY_PERIODIC,50,ctx.physics.dof,2,NULL,&da);
683: DMSetFromOptions(da);
684: DMSetUp(da);
685: /* Inform the DMDA of the field names provided by the physics. */
686: /* The names will be shown in the title bars when run with -ts_monitor_draw_solution */
687: for (i=0; i<ctx.physics.dof; i++) {
688: DMDASetFieldName(da,i,ctx.physics.fieldname[i]);
689: }
690: DMDAGetInfo(da,0,&Mx,0,0,0,0,0,&dof,0,0,0,0,0);
691: DMDAGetCorners(da,&xs,0,0,&xm,0,0);
693: /* Set coordinates of cell centers */
694: DMDASetUniformCoordinates(da,ctx.xmin+0.5*(ctx.xmax-ctx.xmin)/Mx,ctx.xmax+0.5*(ctx.xmax-ctx.xmin)/Mx,0,0,0,0);
696: /* Allocate work space for the Finite Volume solver (so it doesn't have to be reallocated on each function evaluation) */
697: PetscMalloc3(dof,&ctx.u,dof,&ctx.flux,dof,&ctx.speeds);
699: /* Create a vector to store the solution and to save the initial state */
700: DMCreateGlobalVector(da,&X);
701: VecDuplicate(X,&X0);
702: VecDuplicate(X,&R);
704: /* create index for slow parts and fast parts*/
705: count_slow = Mx/(1+ctx.hratio);
707: count_fast = Mx-count_slow;
708: ctx.sf = count_slow/2;
709: ctx.fs = ctx.sf + count_fast;
710: PetscMalloc1(xm*dof,&index_slow);
711: PetscMalloc1(xm*dof,&index_fast);
712: for (i=xs; i<xs+xm; i++) {
713: if (i < count_slow/2 || i > count_slow/2+count_fast-1)
714: for (k=0; k<dof; k++) index_slow[islow++] = i*dof+k;
715: else
716: for (k=0; k<dof; k++) index_fast[ifast++] = i*dof+k;
717: }
718: ISCreateGeneral(PETSC_COMM_WORLD,islow,index_slow,PETSC_COPY_VALUES,&ctx.iss);
719: ISCreateGeneral(PETSC_COMM_WORLD,ifast,index_fast,PETSC_COPY_VALUES,&ctx.isf);
721: /* Create a time-stepping object */
722: TSCreate(comm,&ts);
723: TSSetDM(ts,da);
724: TSSetRHSFunction(ts,R,FVRHSFunction,&ctx);
725: TSRHSSplitSetIS(ts,"slow",ctx.iss);
726: TSRHSSplitSetIS(ts,"fast",ctx.isf);
727: TSRHSSplitSetRHSFunction(ts,"slow",NULL,FVRHSFunctionslow,&ctx);
728: TSRHSSplitSetRHSFunction(ts,"fast",NULL,FVRHSFunctionfast,&ctx);
730: TSSetType(ts,TSMPRK);
731: TSSetMaxTime(ts,10);
732: TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);
734: /* Compute initial conditions and starting time step */
735: FVSample(&ctx,da,0,X0);
736: FVRHSFunction(ts,0,X0,X,(void*)&ctx); /* Initial function evaluation, only used to determine max speed */
737: VecCopy(X0,X); /* The function value was not used so we set X=X0 again */
738: TSSetTimeStep(ts,ctx.cfl/ctx.cfl_idt);
739: TSSetFromOptions(ts); /* Take runtime options */
740: SolutionStatsView(da,X,PETSC_VIEWER_STDOUT_WORLD);
741: {
742: PetscInt steps;
743: PetscScalar mass_initial,mass_final,mass_difference,mass_differenceg;
744: const PetscScalar *ptr_X,*ptr_X0;
745: const PetscReal hs = (ctx.xmax-ctx.xmin)/2.0/count_slow;
746: const PetscReal hf = (ctx.xmax-ctx.xmin)/2.0/count_fast;
747: TSSolve(ts,X);
748: TSGetSolveTime(ts,&ptime);
749: TSGetStepNumber(ts,&steps);
750: /* calculate the total mass at initial time and final time */
751: mass_initial = 0.0;
752: mass_final = 0.0;
753: DMDAVecGetArrayRead(da,X0,(void*)&ptr_X0);
754: DMDAVecGetArrayRead(da,X,(void*)&ptr_X);
755: for (i=xs; i<xs+xm; i++) {
756: if (i < ctx.sf || i > ctx.fs-1) {
757: for (k=0; k<dof; k++) {
758: mass_initial = mass_initial+hs*ptr_X0[i*dof+k];
759: mass_final = mass_final+hs*ptr_X[i*dof+k];
760: }
761: } else {
762: for (k=0; k<dof; k++) {
763: mass_initial = mass_initial+hf*ptr_X0[i*dof+k];
764: mass_final = mass_final+hf*ptr_X[i*dof+k];
765: }
766: }
767: }
768: DMDAVecRestoreArrayRead(da,X0,(void*)&ptr_X0);
769: DMDAVecRestoreArrayRead(da,X,(void*)&ptr_X);
770: mass_difference = mass_final-mass_initial;
771: MPI_Allreduce(&mass_difference,&mass_differenceg,1,MPIU_SCALAR,MPIU_SUM,comm);
772: PetscPrintf(comm,"Mass difference %g\n",(double)mass_differenceg);
773: PetscPrintf(comm,"Final time %g, steps %D\n",(double)ptime,steps);
774: if (ctx.exact) {
775: PetscReal nrm1 = 0;
776: SolutionErrorNorms(&ctx,da,ptime,X,&nrm1);
777: PetscPrintf(comm,"Error ||x-x_e||_1 %g\n",(double)nrm1);
778: }
779: if (ctx.simulation) {
780: PetscReal nrm1 = 0;
781: PetscViewer fd;
782: char filename[PETSC_MAX_PATH_LEN] = "binaryoutput";
783: Vec XR;
784: PetscBool flg;
785: const PetscScalar *ptr_XR;
786: PetscOptionsGetString(NULL,NULL,"-f",filename,sizeof(filename),&flg);
788: PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_READ,&fd);
789: VecDuplicate(X0,&XR);
790: VecLoad(XR,fd);
791: PetscViewerDestroy(&fd);
792: VecGetArrayRead(X,&ptr_X);
793: VecGetArrayRead(XR,&ptr_XR);
794: for (i=0; i<Mx; i++) {
795: if (i < count_slow/2 || i > count_slow/2+count_fast-1) nrm1 = nrm1 + hs*PetscAbs(ptr_X[i]-ptr_XR[i]);
796: else nrm1 = nrm1 + hf*PetscAbs(ptr_X[i]-ptr_XR[i]);
797: }
798: VecRestoreArrayRead(X,&ptr_X);
799: VecRestoreArrayRead(XR,&ptr_XR);
800: PetscPrintf(comm,"Error ||x-x_e||_1 %g\n",(double)nrm1);
801: VecDestroy(&XR);
802: }
803: }
805: SolutionStatsView(da,X,PETSC_VIEWER_STDOUT_WORLD);
806: if (draw & 0x1) VecView(X0,PETSC_VIEWER_DRAW_WORLD);
807: if (draw & 0x2) VecView(X,PETSC_VIEWER_DRAW_WORLD);
808: if (draw & 0x4) {
809: Vec Y;
810: VecDuplicate(X,&Y);
811: FVSample(&ctx,da,ptime,Y);
812: VecAYPX(Y,-1,X);
813: VecView(Y,PETSC_VIEWER_DRAW_WORLD);
814: VecDestroy(&Y);
815: }
817: if (view_final) {
818: PetscViewer viewer;
819: PetscViewerASCIIOpen(PETSC_COMM_WORLD,final_fname,&viewer);
820: PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_MATLAB);
821: VecView(X,viewer);
822: PetscViewerPopFormat(viewer);
823: PetscViewerDestroy(&viewer);
824: }
826: /* Clean up */
827: (*ctx.physics.destroy)(ctx.physics.user);
828: for (i=0; i<ctx.physics.dof; i++) PetscFree(ctx.physics.fieldname[i]);
829: PetscFree3(ctx.u,ctx.flux,ctx.speeds);
830: ISDestroy(&ctx.iss);
831: ISDestroy(&ctx.isf);
832: VecDestroy(&X);
833: VecDestroy(&X0);
834: VecDestroy(&R);
835: DMDestroy(&da);
836: TSDestroy(&ts);
837: PetscFree(index_slow);
838: PetscFree(index_fast);
839: PetscFunctionListDestroy(&physics);
840: PetscFinalize();
841: return 0;
842: }
844: /*TEST
846: build:
847: requires: !complex
849: test:
850: 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
852: test:
853: suffix: 2
854: 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
855: output_file: output/ex7_1.out
857: test:
858: suffix: 3
859: 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
861: test:
862: suffix: 4
863: 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
864: output_file: output/ex7_3.out
866: test:
867: suffix: 5
868: nsize: 2
869: 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
870: output_file: output/ex7_3.out
871: TEST*/