Actual source code: ex7.c
petsc-3.12.0 2019-09-29
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,1,"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,1,"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 to 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: if (1) {
289: PetscPrintf(ctx->comm,"Stability constraint exceeded at t=%g, dt %g > %g\n",(double)tnow,(double)dt,(double)(0.5/ctx->cfl_idt));
290: } else SETERRQ2(PETSC_COMM_SELF,1,"Stability constraint exceeded, %g > %g",(double)dt,(double)(ctx->cfl/ctx->cfl_idt));
291: }
292: }
293: PetscFree4(r,min,alpha,gamma);
294: return(0);
295: }
297: static PetscErrorCode FVRHSFunctionslow(TS ts,PetscReal time,Vec X,Vec F,void *vctx)
298: {
299: FVCtx *ctx = (FVCtx*)vctx;
300: PetscErrorCode ierr;
301: PetscInt i,j,Mx,dof,xs,xm,islow = 0,sf = ctx->sf,fs = ctx->fs;
302: PetscReal hf,hs;
303: PetscScalar *x,*f,*r,*min,*alpha,*gamma;
304: Vec Xloc;
305: DM da;
308: TSGetDM(ts,&da);
309: DMGetLocalVector(da,&Xloc); /* Xloc contains ghost points */
310: DMDAGetInfo(da,0,&Mx,0,0,0,0,0,&dof,0,0,0,0,0); /* Mx is the number of center points */
311: hs = (ctx->xmax-ctx->xmin)/2.0*(ctx->hratio+1.0)/Mx;
312: hf = (ctx->xmax-ctx->xmin)/2.0*(1.0+1.0/ctx->hratio)/Mx;
313: DMGlobalToLocalBegin(da,X,INSERT_VALUES,Xloc); /* X is solution vector which does not contain ghost points */
314: DMGlobalToLocalEnd (da,X,INSERT_VALUES,Xloc);
315: VecZeroEntries(F); /* F is the right hand side function corresponds to center points */
316: DMDAVecGetArray(da,Xloc,&x);
317: VecGetArray(F,&f);
318: DMDAGetCorners(da,&xs,0,0,&xm,0,0);
319: PetscMalloc4(dof,&r,dof,&min,dof,&alpha,dof,&gamma);
321: if (ctx->bctype == FVBC_OUTFLOW) {
322: for (i=xs-2; i<0; i++) {
323: for (j=0; j<dof; j++) x[i*dof+j] = x[j];
324: }
325: for (i=Mx; i<xs+xm+2; i++) {
326: for (j=0; j<dof; j++) x[i*dof+j] = x[(xs+xm-1)*dof+j];
327: }
328: }
330: for (i=xs; i<xs+xm+1; i++) {
331: PetscReal maxspeed;
332: PetscScalar *u;
333: if (i < sf) {
334: u = &ctx->u[0];
335: alpha[0] = 1.0/6.0;
336: gamma[0] = 1.0/3.0;
337: for (j=0; j<dof; j++) {
338: r[j] = (x[i*dof+j] - x[(i-1)*dof+j])/(x[(i-1)*dof+j]-x[(i-2)*dof+j]);
339: min[j] = PetscMin(r[j],2.0);
340: 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]);
341: }
342: (*ctx->physics.flux)(ctx->physics.user,u,ctx->flux,&maxspeed);
343: if (i > xs) {
344: for (j=0; j<dof; j++) f[(i-1)*dof+j] -= ctx->flux[j]/hs;
345: }
346: if (i < xs+xm) {
347: for (j=0; j<dof; j++) f[i*dof+j] += ctx->flux[j]/hs;
348: islow++;
349: }
350: } else if (i == sf) {
351: u = &ctx->u[0];
352: alpha[0] = hs*hf/(hs+hs)/(hs+hs+hf);
353: gamma[0] = hs*(hs+hs)/(hs+hf)/(hs+hs+hf);
354: for (j=0; j<dof; j++) {
355: r[j] = (x[i*dof+j] - x[(i-1)*dof+j])/(x[(i-1)*dof+j]-x[(i-2)*dof+j]);
356: min[j] = PetscMin(r[j],2.0);
357: 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]);
358: }
359: (*ctx->physics.flux)(ctx->physics.user,u,ctx->flux,&maxspeed);
360: if (i < xs+xm) {
361: for (j=0; j<dof; j++) f[(i-1)*dof+j] -= ctx->flux[j]/hs;
362: islow++;
363: }
364: } else if (i == fs) {
365: u = &ctx->u[0];
366: alpha[0] = hf*hs/(hf+hf)/(hf+hf+hs);
367: gamma[0] = hf*(hf+hf)/(hf+hs)/(hf+hf+hs);
368: for (j=0; j<dof; j++) {
369: r[j] = (x[i*dof+j] - x[(i-1)*dof+j])/(x[(i-1)*dof+j]-x[(i-2)*dof+j]);
370: min[j] = PetscMin(r[j],2.0);
371: 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]);
372: }
373: (*ctx->physics.flux)(ctx->physics.user,u,ctx->flux,&maxspeed);
374: if (i < xs+xm) {
375: for (j=0; j<dof; j++) f[i*dof+j-fs+sf] += ctx->flux[j]/hs;
376: islow++;
377: }
378: } else if (i == fs+1) {
379: u = &ctx->u[0];
380: alpha[0] = hs*hs/(hf+hs)/(hf+hs+hs);
381: gamma[0] = hs*(hf+hs)/(hs+hs)/(hf+hs+hs);
382: for (j=0; j<dof; j++) {
383: r[j] = (x[i*dof+j] - x[(i-1)*dof+j])/(x[(i-1)*dof+j]-x[(i-2)*dof+j]);
384: min[j] = PetscMin(r[j],2.0);
385: 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]);
386: }
387: (*ctx->physics.flux)(ctx->physics.user,u,ctx->flux,&maxspeed);
388: if (i > xs) {
389: for (j=0; j<dof; j++) f[(i-fs+sf-1)*dof+j] -= ctx->flux[j]/hs;
390: }
391: if (i < xs+xm) {
392: for (j=0; j<dof; j++) f[(i-fs+sf)*dof+j] += ctx->flux[j]/hs;
393: }
394: } else if (i > fs+1) {
395: u = &ctx->u[0];
396: alpha[0] = 1.0/6.0;
397: gamma[0] = 1.0/3.0;
398: for (j=0; j<dof; j++) {
399: r[j] = (x[i*dof+j] - x[(i-1)*dof+j])/(x[(i-1)*dof+j]-x[(i-2)*dof+j]);
400: min[j] = PetscMin(r[j],2.0);
401: 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]);
402: }
403: (*ctx->physics.flux)(ctx->physics.user,u,ctx->flux,&maxspeed);
404: if (i > xs) {
405: for (j=0; j<dof; j++) f[(i-fs+sf-1)*dof+j] -= ctx->flux[j]/hs;
406: }
407: if (i < xs+xm) {
408: for (j=0; j<dof; j++) f[(i-fs+sf)*dof+j] += ctx->flux[j]/hs;
409: }
410: }
411: }
412: DMDAVecRestoreArray(da,Xloc,&x);
413: VecRestoreArray(F,&f);
414: DMRestoreLocalVector(da,&Xloc);
415: PetscFree4(r,min,alpha,gamma);
416: return(0);
417: }
419: static PetscErrorCode FVRHSFunctionfast(TS ts,PetscReal time,Vec X,Vec F,void *vctx)
420: {
421: FVCtx *ctx = (FVCtx*)vctx;
423: PetscInt i,j,Mx,dof,xs,xm,ifast = 0,sf = ctx->sf,fs = ctx->fs;
424: PetscReal hf,hs;
425: PetscScalar *x,*f,*r,*min,*alpha,*gamma;
426: Vec Xloc;
427: DM da;
430: TSGetDM(ts,&da);
431: DMGetLocalVector(da,&Xloc); /* Xloc contains ghost points */
432: DMDAGetInfo(da,0,&Mx,0,0,0,0,0,&dof,0,0,0,0,0); /* Mx is the number of center points */
433: hs = (ctx->xmax-ctx->xmin)/2.0*(ctx->hratio+1.0)/Mx;
434: hf = (ctx->xmax-ctx->xmin)/2.0*(1.0+1.0/ctx->hratio)/Mx;
435: DMGlobalToLocalBegin(da,X,INSERT_VALUES,Xloc); /* X is solution vector which does not contain ghost points */
436: DMGlobalToLocalEnd(da,X,INSERT_VALUES,Xloc);
437: VecZeroEntries(F); /* F is the right hand side function corresponds to center points */
438: DMDAVecGetArray(da,Xloc,&x);
439: VecGetArray(F,&f);
440: DMDAGetCorners(da,&xs,0,0,&xm,0,0);
441: PetscMalloc4(dof,&r,dof,&min,dof,&alpha,dof,&gamma);
443: if (ctx->bctype == FVBC_OUTFLOW) {
444: for (i=xs-2; i<0; i++) {
445: for (j=0; j<dof; j++) x[i*dof+j] = x[j];
446: }
447: for (i=Mx; i<xs+xm+2; i++) {
448: for (j=0; j<dof; j++) x[i*dof+j] = x[(xs+xm-1)*dof+j];
449: }
450: }
452: for (i=xs; i<xs+xm+1; i++) {
453: PetscReal maxspeed;
454: PetscScalar *u;
455: if (i == sf) {
456: u = &ctx->u[0];
457: alpha[0] = hs*hf/(hs+hs)/(hs+hs+hf);
458: gamma[0] = hs*(hs+hs)/(hs+hf)/(hs+hs+hf);
459: for (j=0; j<dof; j++) {
460: r[j] = (x[i*dof+j]-x[(i-1)*dof+j])/(x[(i-1)*dof+j]-x[(i-2)*dof+j]);
461: min[j] = PetscMin(r[j],2.0);
462: 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]);
463: }
464: (*ctx->physics.flux)(ctx->physics.user,u,ctx->flux,&maxspeed);
465: if (i < xs+xm) {
466: for (j=0; j<dof; j++) f[(i-sf)*dof+j] += ctx->flux[j]/hf;
467: ifast++;
468: }
469: } else if (i == sf+1) {
470: u = &ctx->u[0];
471: alpha[0] = hf*hf/(hs+hf)/(hs+hf+hf);
472: gamma[0] = hf*(hs+hf)/(hf+hf)/(hs+hf+hf);
473: for (j=0; j<dof; j++) {
474: r[j] = (x[i*dof+j]-x[(i-1)*dof+j])/(x[(i-1)*dof+j]-x[(i-2)*dof+j]);
475: min[j] = PetscMin(r[j],2.0);
476: 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]);
477: }
478: (*ctx->physics.flux)(ctx->physics.user,u,ctx->flux,&maxspeed);
479: if (i > xs) {
480: for (j=0; j<dof; j++) f[(i-sf-1)*dof+j] -= ctx->flux[j]/hf;
481: }
482: if (i < xs+xm) {
483: for (j=0; j<dof; j++) f[(i-sf)*dof+j] += ctx->flux[j]/hf;
484: ifast++;
485: }
486: } else if (i > sf+1 && i < fs) {
487: u = &ctx->u[0];
488: alpha[0] = 1.0/6.0;
489: gamma[0] = 1.0/3.0;
490: for (j=0; j<dof; j++) {
491: r[j] = (x[i*dof+j]-x[(i-1)*dof+j])/(x[(i-1)*dof+j]-x[(i-2)*dof+j]);
492: min[j] = PetscMin(r[j],2.0);
493: 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]);
494: }
495: (*ctx->physics.flux)(ctx->physics.user,u,ctx->flux,&maxspeed);
496: if (i > xs) {
497: for (j=0; j<dof; j++) f[(i-sf-1)*dof+j] -= ctx->flux[j]/hf;
498: }
499: if (i < xs+xm) {
500: for (j=0; j<dof; j++) f[(i-sf)*dof+j] += ctx->flux[j]/hf;
501: ifast++;
502: }
503: } else if (i == fs) {
504: u = &ctx->u[0];
505: alpha[0] = hf*hs/(hf+hf)/(hf+hf+hs);
506: gamma[0] = hf*(hf+hf)/(hf+hs)/(hf+hf+hs);
507: for (j=0; j<dof; j++) {
508: r[j] = (x[i*dof+j]-x[(i-1)*dof+j])/(x[(i-1)*dof+j]-x[(i-2)*dof+j]);
509: min[j] = PetscMin(r[j],2.0);
510: 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]);
511: }
512: (*ctx->physics.flux)(ctx->physics.user,u,ctx->flux,&maxspeed);
513: if (i>xs) {
514: for (j=0; j<dof; j++) f[(i-sf-1)*dof+j] -= ctx->flux[j]/hf;
515: }
516: }
517: }
518: DMDAVecRestoreArray(da,Xloc,&x);
519: VecRestoreArray(F,&f);
520: DMRestoreLocalVector(da,&Xloc);
521: PetscFree4(r,min,alpha,gamma);
522: return(0);
523: }
525: /* --------------------------------- Finite Volume Solver for slow components ----------------------------------- */
527: PetscErrorCode FVSample(FVCtx *ctx,DM da,PetscReal time,Vec U)
528: {
530: PetscScalar *u,*uj,xj,xi;
531: PetscInt i,j,k,dof,xs,xm,Mx,count_slow,count_fast;
532: const PetscInt N=200;
535: if (!ctx->physics.sample) SETERRQ(PETSC_COMM_SELF,1,"Physics has not provided a sampling function");
536: DMDAGetInfo(da,0,&Mx,0,0,0,0,0,&dof,0,0,0,0,0);
537: DMDAGetCorners(da,&xs,0,0,&xm,0,0);
538: DMDAVecGetArray(da,U,&u);
539: PetscMalloc1(dof,&uj);
540: const PetscReal hs = (ctx->xmax-ctx->xmin)/2.0*(ctx->hratio+1.0)/Mx;
541: const PetscReal hf = (ctx->xmax-ctx->xmin)/2.0*(1.0+1.0/ctx->hratio)/Mx;
542: count_slow = Mx/(1+ctx->hratio);
543: count_fast = Mx-count_slow;
544: for (i=xs; i<xs+xm; i++) {
545: if (i*hs+0.5*hs<(ctx->xmax-ctx->xmin)*0.25) {
546: xi = ctx->xmin+0.5*hs+i*hs;
547: /* Integrate over cell i using trapezoid rule with N points. */
548: for (k=0; k<dof; k++) u[i*dof+k] = 0;
549: for (j=0; j<N+1; j++) {
550: xj = xi+hs*(j-N/2)/(PetscReal)N;
551: (*ctx->physics.sample)(ctx->physics.user,ctx->initial,ctx->bctype,ctx->xmin,ctx->xmax,time,xj,uj);
552: for (k=0; k<dof; k++) u[i*dof+k] += ((j==0 || j==N) ? 0.5 : 1.0)*uj[k]/N;
553: }
554: } else if ((ctx->xmax-ctx->xmin)*0.25+(i-count_slow/2)*hf+0.5*hf<(ctx->xmax-ctx->xmin)*0.75) {
555: xi = ctx->xmin+(ctx->xmax-ctx->xmin)*0.25+0.5*hf+(i-count_slow/2)*hf;
556: /* Integrate over cell i using trapezoid rule with N points. */
557: for (k=0; k<dof; k++) u[i*dof+k] = 0;
558: for (j=0; j<N+1; j++) {
559: xj = xi+hf*(j-N/2)/(PetscReal)N;
560: (*ctx->physics.sample)(ctx->physics.user,ctx->initial,ctx->bctype,ctx->xmin,ctx->xmax,time,xj,uj);
561: for (k=0; k<dof; k++) u[i*dof+k] += ((j==0 || j==N) ? 0.5 : 1.0)*uj[k]/N;
562: }
563: } else {
564: xi = ctx->xmin+(ctx->xmax-ctx->xmin)*0.75+0.5*hs+(i-count_slow/2-count_fast)*hs;
565: /* Integrate over cell i using trapezoid rule with N points. */
566: for (k=0; k<dof; k++) u[i*dof+k] = 0;
567: for (j=0; j<N+1; j++) {
568: xj = xi+hs*(j-N/2)/(PetscReal)N;
569: (*ctx->physics.sample)(ctx->physics.user,ctx->initial,ctx->bctype,ctx->xmin,ctx->xmax,time,xj,uj);
570: for (k=0; k<dof; k++) u[i*dof+k] += ((j==0 || j==N) ? 0.5 : 1.0)*uj[k]/N;
571: }
572: }
573: }
574: DMDAVecRestoreArray(da,U,&u);
575: PetscFree(uj);
576: return(0);
577: }
579: static PetscErrorCode SolutionStatsView(DM da,Vec X,PetscViewer viewer)
580: {
581: PetscErrorCode ierr;
582: PetscReal xmin,xmax;
583: PetscScalar sum,tvsum,tvgsum;
584: const PetscScalar *x;
585: PetscInt imin,imax,Mx,i,j,xs,xm,dof;
586: Vec Xloc;
587: PetscBool iascii;
590: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
591: if (iascii) {
592: /* PETSc lacks a function to compute total variation norm (difficult in multiple dimensions), we do it here */
593: DMGetLocalVector(da,&Xloc);
594: DMGlobalToLocalBegin(da,X,INSERT_VALUES,Xloc);
595: DMGlobalToLocalEnd(da,X,INSERT_VALUES,Xloc);
596: DMDAVecGetArrayRead(da,Xloc,(void*)&x);
597: DMDAGetCorners(da,&xs,0,0,&xm,0,0);
598: DMDAGetInfo(da,0,&Mx,0,0,0,0,0,&dof,0,0,0,0,0);
599: tvsum = 0;
600: for (i=xs; i<xs+xm; i++) {
601: for (j=0; j<dof; j++) tvsum += PetscAbsScalar(x[i*dof+j]-x[(i-1)*dof+j]);
602: }
603: MPI_Allreduce(&tvsum,&tvgsum,1,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)da));
604: DMDAVecRestoreArrayRead(da,Xloc,(void*)&x);
605: DMRestoreLocalVector(da,&Xloc);
607: VecMin(X,&imin,&xmin);
608: VecMax(X,&imax,&xmax);
609: VecSum(X,&sum);
610: 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));
611: } else SETERRQ(PETSC_COMM_SELF,1,"Viewer type not supported");
612: return(0);
613: }
615: static PetscErrorCode SolutionErrorNorms(FVCtx *ctx,DM da,PetscReal t,Vec X,PetscReal *nrm1)
616: {
617: PetscErrorCode ierr;
618: Vec Y;
619: PetscInt i,Mx,count_slow=0,count_fast=0;
620: const PetscScalar *ptr_X,*ptr_Y;
623: VecGetSize(X,&Mx);
624: VecDuplicate(X,&Y);
625: FVSample(ctx,da,t,Y);
626: const PetscReal hs = (ctx->xmax-ctx->xmin)/2.0*(ctx->hratio+1.0)/Mx;
627: const PetscReal hf = (ctx->xmax-ctx->xmin)/2.0*(1.0+1.0/ctx->hratio)/Mx;
628: count_slow = (PetscReal)Mx/(1.0+ctx->hratio);
629: count_fast = Mx-count_slow;
630: VecGetArrayRead(X,&ptr_X);
631: VecGetArrayRead(Y,&ptr_Y);
632: for (i=0; i<Mx; i++) {
633: if (i < count_slow/2 || i > count_slow/2+count_fast-1) *nrm1 += hs*PetscAbs(ptr_X[i]-ptr_Y[i]);
634: else *nrm1 += hf*PetscAbs(ptr_X[i]-ptr_Y[i]);
635: }
636: VecRestoreArrayRead(X,&ptr_X);
637: VecRestoreArrayRead(Y,&ptr_Y);
638: VecDestroy(&Y);
639: return(0);
640: }
642: int main(int argc,char *argv[])
643: {
644: char physname[256] = "advect",final_fname[256] = "solution.m";
645: PetscFunctionList physics = 0;
646: MPI_Comm comm;
647: TS ts;
648: DM da;
649: Vec X,X0,R;
650: FVCtx ctx;
651: PetscInt i,k,dof,xs,xm,Mx,draw = 0,count_slow,count_fast,islow = 0,ifast = 0,*index_slow,*index_fast;
652: PetscBool view_final = PETSC_FALSE;
653: PetscReal ptime;
654: PetscErrorCode ierr;
656: PetscInitialize(&argc,&argv,0,help);if (ierr) return ierr;
657: comm = PETSC_COMM_WORLD;
658: PetscMemzero(&ctx,sizeof(ctx));
660: /* Register physical models to be available on the command line */
661: PetscFunctionListAdd(&physics,"advect",PhysicsCreate_Advect);
663: ctx.comm = comm;
664: ctx.cfl = 0.9;
665: ctx.bctype = FVBC_PERIODIC;
666: ctx.xmin = -1.0;
667: ctx.xmax = 1.0;
668: PetscOptionsBegin(comm,NULL,"Finite Volume solver options","");
669: PetscOptionsReal("-xmin","X min","",ctx.xmin,&ctx.xmin,NULL);
670: PetscOptionsReal("-xmax","X max","",ctx.xmax,&ctx.xmax,NULL);
671: PetscOptionsInt("-draw","Draw solution vector, bitwise OR of (1=initial,2=final,4=final error)","",draw,&draw,NULL);
672: PetscOptionsString("-view_final","Write final solution in ASCII MATLAB format to given file name","",final_fname,final_fname,sizeof(final_fname),&view_final);
673: PetscOptionsInt("-initial","Initial condition (depends on the physics)","",ctx.initial,&ctx.initial,NULL);
674: PetscOptionsBool("-exact","Compare errors with exact solution","",ctx.exact,&ctx.exact,NULL);
675: PetscOptionsBool("-simulation","Compare errors with reference solution","",ctx.simulation,&ctx.simulation,NULL);
676: PetscOptionsReal("-cfl","CFL number to time step at","",ctx.cfl,&ctx.cfl,NULL);
677: PetscOptionsEnum("-bc_type","Boundary condition","",FVBCTypes,(PetscEnum)ctx.bctype,(PetscEnum*)&ctx.bctype,NULL);
678: PetscOptionsInt("-hratio","Spacing ratio","",ctx.hratio,&ctx.hratio,NULL);
679: PetscOptionsEnd();
681: /* Choose the physics from the list of registered models */
682: {
683: PetscErrorCode (*r)(FVCtx*);
684: PetscFunctionListFind(physics,physname,&r);
685: if (!r) SETERRQ1(PETSC_COMM_SELF,1,"Physics '%s' not found",physname);
686: /* Create the physics, will set the number of fields and their names */
687: (*r)(&ctx);
688: }
690: /* Create a DMDA to manage the parallel grid */
691: DMDACreate1d(comm,DM_BOUNDARY_PERIODIC,50,ctx.physics.dof,2,NULL,&da);
692: DMSetFromOptions(da);
693: DMSetUp(da);
694: /* Inform the DMDA of the field names provided by the physics. */
695: /* The names will be shown in the title bars when run with -ts_monitor_draw_solution */
696: for (i=0; i<ctx.physics.dof; i++) {
697: DMDASetFieldName(da,i,ctx.physics.fieldname[i]);
698: }
699: DMDAGetInfo(da,0,&Mx,0,0,0,0,0,&dof,0,0,0,0,0);
700: DMDAGetCorners(da,&xs,0,0,&xm,0,0);
702: /* Set coordinates of cell centers */
703: DMDASetUniformCoordinates(da,ctx.xmin+0.5*(ctx.xmax-ctx.xmin)/Mx,ctx.xmax+0.5*(ctx.xmax-ctx.xmin)/Mx,0,0,0,0);
705: /* Allocate work space for the Finite Volume solver (so it doesn't have to be reallocated on each function evaluation) */
706: PetscMalloc3(dof,&ctx.u,dof,&ctx.flux,dof,&ctx.speeds);
708: /* Create a vector to store the solution and to save the initial state */
709: DMCreateGlobalVector(da,&X);
710: VecDuplicate(X,&X0);
711: VecDuplicate(X,&R);
713: /* create index for slow parts and fast parts*/
714: count_slow = Mx/(1+ctx.hratio);
715: if (count_slow%2) SETERRQ(PETSC_COMM_WORLD,1,"Please adjust grid size Mx (-da_grid_x) and hratio (-hratio) so that Mx/(1+hartio) is even");
716: count_fast = Mx-count_slow;
717: ctx.sf = count_slow/2;
718: ctx.fs = ctx.sf + count_fast;
719: PetscMalloc1(xm*dof,&index_slow);
720: PetscMalloc1(xm*dof,&index_fast);
721: for (i=xs; i<xs+xm; i++) {
722: if (i < count_slow/2 || i > count_slow/2+count_fast-1)
723: for (k=0; k<dof; k++) index_slow[islow++] = i*dof+k;
724: else
725: for (k=0; k<dof; k++) index_fast[ifast++] = i*dof+k;
726: }
727: ISCreateGeneral(PETSC_COMM_WORLD,islow,index_slow,PETSC_COPY_VALUES,&ctx.iss);
728: ISCreateGeneral(PETSC_COMM_WORLD,ifast,index_fast,PETSC_COPY_VALUES,&ctx.isf);
730: /* Create a time-stepping object */
731: TSCreate(comm,&ts);
732: TSSetDM(ts,da);
733: TSSetRHSFunction(ts,R,FVRHSFunction,&ctx);
734: TSRHSSplitSetIS(ts,"slow",ctx.iss);
735: TSRHSSplitSetIS(ts,"fast",ctx.isf);
736: TSRHSSplitSetRHSFunction(ts,"slow",NULL,FVRHSFunctionslow,&ctx);
737: TSRHSSplitSetRHSFunction(ts,"fast",NULL,FVRHSFunctionfast,&ctx);
739: TSSetType(ts,TSMPRK);
740: TSSetMaxTime(ts,10);
741: TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);
743: /* Compute initial conditions and starting time step */
744: FVSample(&ctx,da,0,X0);
745: FVRHSFunction(ts,0,X0,X,(void*)&ctx); /* Initial function evaluation, only used to determine max speed */
746: VecCopy(X0,X); /* The function value was not used so we set X=X0 again */
747: TSSetTimeStep(ts,ctx.cfl/ctx.cfl_idt);
748: TSSetFromOptions(ts); /* Take runtime options */
749: SolutionStatsView(da,X,PETSC_VIEWER_STDOUT_WORLD);
750: {
751: PetscInt steps;
752: PetscScalar mass_initial,mass_final,mass_difference,mass_differenceg;
753: const PetscScalar *ptr_X,*ptr_X0;
754: const PetscReal hs = (ctx.xmax-ctx.xmin)/2.0*(ctx.hratio+1.0)/Mx;
755: const PetscReal hf = (ctx.xmax-ctx.xmin)/2.0*(1.0+1.0/ctx.hratio)/Mx;
756: TSSolve(ts,X);
757: TSGetSolveTime(ts,&ptime);
758: TSGetStepNumber(ts,&steps);
759: /* calculate the total mass at initial time and final time */
760: mass_initial = 0.0;
761: mass_final = 0.0;
762: DMDAVecGetArrayRead(da,X0,(void*)&ptr_X0);
763: DMDAVecGetArrayRead(da,X,(void*)&ptr_X);
764: for(i=0; i<Mx; i++) {
765: if(i < count_slow/2 || i > count_slow/2+count_fast-1){
766: mass_initial = mass_initial+hs*ptr_X0[i];
767: mass_final = mass_final+hs*ptr_X[i];
768: } else {
769: mass_initial = mass_initial+hf*ptr_X0[i];
770: mass_final = mass_final+hf*ptr_X[i];
771: }
772: }
773: DMDAVecRestoreArrayRead(da,X0,(void*)&ptr_X0);
774: DMDAVecRestoreArrayRead(da,X,(void*)&ptr_X);
775: mass_difference = mass_final-mass_initial;
776: MPI_Allreduce(&mass_difference,&mass_differenceg,1,MPIU_SCALAR,MPIU_SUM,comm);
777: PetscPrintf(comm,"Mass difference %g\n",(double)mass_differenceg);
778: PetscPrintf(comm,"Final time %g, steps %D\n",(double)ptime,steps);
779: if (ctx.exact) {
780: PetscReal nrm1=0;
781: SolutionErrorNorms(&ctx,da,ptime,X,&nrm1);
782: PetscPrintf(comm,"Error ||x-x_e||_1 %g\n",(double)nrm1);
783: }
784: if (ctx.simulation) {
785: PetscReal nrm1=0;
786: PetscViewer fd;
787: char filename[PETSC_MAX_PATH_LEN] = "binaryoutput";
788: Vec XR;
789: PetscBool flg;
790: const PetscScalar *ptr_XR;
791: PetscOptionsGetString(NULL,NULL,"-f",filename,sizeof(filename),&flg);
792: if (!flg) SETERRQ(PETSC_COMM_WORLD,1,"Must indicate binary file with the -f option");
793: PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_READ,&fd);
794: VecDuplicate(X0,&XR);
795: VecLoad(XR,fd);
796: PetscViewerDestroy(&fd);
797: VecGetArrayRead(X,&ptr_X);
798: VecGetArrayRead(XR,&ptr_XR);
799: for(i=0; i<Mx; i++) {
800: if(i < count_slow/2 || i > count_slow/2+count_fast-1) nrm1 = nrm1 + hs*PetscAbs(ptr_X[i]-ptr_XR[i]);
801: else nrm1 = nrm1 + hf*PetscAbs(ptr_X[i]-ptr_XR[i]);
802: }
803: VecRestoreArrayRead(X,&ptr_X);
804: VecRestoreArrayRead(XR,&ptr_XR);
805: PetscPrintf(comm,"Error ||x-x_e||_1 %g\n",(double)nrm1);
806: VecDestroy(&XR);
807: }
808: }
810: SolutionStatsView(da,X,PETSC_VIEWER_STDOUT_WORLD);
811: if (draw & 0x1) { VecView(X0,PETSC_VIEWER_DRAW_WORLD); }
812: if (draw & 0x2) { VecView(X,PETSC_VIEWER_DRAW_WORLD); }
813: if (draw & 0x4) {
814: Vec Y;
815: VecDuplicate(X,&Y);
816: FVSample(&ctx,da,ptime,Y);
817: VecAYPX(Y,-1,X);
818: VecView(Y,PETSC_VIEWER_DRAW_WORLD);
819: VecDestroy(&Y);
820: }
822: if (view_final) {
823: PetscViewer viewer;
824: PetscViewerASCIIOpen(PETSC_COMM_WORLD,final_fname,&viewer);
825: PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_MATLAB);
826: VecView(X,viewer);
827: PetscViewerPopFormat(viewer);
828: PetscViewerDestroy(&viewer);
829: }
831: /* Clean up */
832: (*ctx.physics.destroy)(ctx.physics.user);
833: for (i=0; i<ctx.physics.dof; i++) {PetscFree(ctx.physics.fieldname[i]);}
834: PetscFree3(ctx.u,ctx.flux,ctx.speeds);
835: ISDestroy(&ctx.iss);
836: ISDestroy(&ctx.isf);
837: VecDestroy(&X);
838: VecDestroy(&X0);
839: VecDestroy(&R);
840: DMDestroy(&da);
841: TSDestroy(&ts);
842: PetscFree(index_slow);
843: PetscFree(index_fast);
844: PetscFunctionListDestroy(&physics);
845: PetscFinalize();
846: return ierr;
847: }
849: /*TEST
851: build:
852: requires: !complex c99
854: test:
855: 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
857: test:
858: suffix: 2
859: 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
860: output_file: output/ex7_1.out
862: test:
863: suffix: 3
864: 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
866: test:
867: suffix: 4
868: 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
869: output_file: output/ex7_3.out
870: TEST*/