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