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