Actual source code: ex5.c
petsc-3.12.0 2019-09-29
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> /* For the Kernel_*_gets_* stuff for BAIJ */
32: #include "finitevolume1d.h"
34: PETSC_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,1,"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,1,"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,1,"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,1,"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;
144: PetscInt i,j,k,Mx,dof,xs,xm,len_slow;
145: PetscReal hx,cfl_idt = 0;
146: PetscScalar *x,*f,*slope;
147: Vec Xloc;
148: DM da;
151: TSGetDM(ts,&da);
152: DMGetLocalVector(da,&Xloc);
153: DMDAGetInfo(da,0, &Mx,0,0, 0,0,0, &dof,0,0,0,0,0);
154: hx = (ctx->xmax-ctx->xmin)/Mx;
155: DMGlobalToLocalBegin(da,X,INSERT_VALUES,Xloc);
156: DMGlobalToLocalEnd(da,X,INSERT_VALUES,Xloc);
157: VecZeroEntries(F);
158: DMDAVecGetArray(da,Xloc,&x);
159: VecGetArray(F,&f);
160: DMDAGetArray(da,PETSC_TRUE,&slope);
161: DMDAGetCorners(da,&xs,0,0,&xm,0,0);
162: ISGetSize(ctx->iss,&len_slow);
164: if (ctx->bctype == FVBC_OUTFLOW) {
165: for (i=xs-2; i<0; i++) {
166: for (j=0; j<dof; j++) x[i*dof+j] = x[j];
167: }
168: for (i=Mx; i<xs+xm+2; i++) {
169: for (j=0; j<dof; j++) x[i*dof+j] = x[(xs+xm-1)*dof+j];
170: }
171: }
172: for (i=xs-1; i<xs+xm+1; i++) {
173: struct _LimitInfo info;
174: PetscScalar *cjmpL,*cjmpR;
175: if (i < len_slow+1) {
176: /* Determine the right eigenvectors R, where A = R \Lambda R^{-1} */
177: (*ctx->physics.characteristic)(ctx->physics.user,dof,&x[i*dof],ctx->R,ctx->Rinv,ctx->speeds,ctx->xmin+hx*i);
178: /* Evaluate jumps across interfaces (i-1, i) and (i, i+1), put in characteristic basis */
179: PetscArrayzero(ctx->cjmpLR,2*dof);
180: cjmpL = &ctx->cjmpLR[0];
181: cjmpR = &ctx->cjmpLR[dof];
182: for (j=0; j<dof; j++) {
183: PetscScalar jmpL,jmpR;
184: jmpL = x[(i+0)*dof+j]-x[(i-1)*dof+j];
185: jmpR = x[(i+1)*dof+j]-x[(i+0)*dof+j];
186: for (k=0; k<dof; k++) {
187: cjmpL[k] += ctx->Rinv[k+j*dof]*jmpL;
188: cjmpR[k] += ctx->Rinv[k+j*dof]*jmpR;
189: }
190: }
191: /* Apply limiter to the left and right characteristic jumps */
192: info.m = dof;
193: info.hx = hx;
194: (*ctx->limit)(&info,cjmpL,cjmpR,ctx->cslope);
195: for (j=0; j<dof; j++) ctx->cslope[j] /= hx; /* rescale to a slope */
196: for (j=0; j<dof; j++) {
197: PetscScalar tmp = 0;
198: for (k=0; k<dof; k++) tmp += ctx->R[j+k*dof]*ctx->cslope[k];
199: slope[i*dof+j] = tmp;
200: }
201: }
202: }
204: for (i=xs; i<xs+xm+1; i++) {
205: PetscReal maxspeed;
206: PetscScalar *uL,*uR;
207: if (i < len_slow) { /* slow parts can be changed based on a */
208: uL = &ctx->uLR[0];
209: uR = &ctx->uLR[dof];
210: for (j=0; j<dof; j++) {
211: uL[j] = x[(i-1)*dof+j]+slope[(i-1)*dof+j]*hx/2;
212: uR[j] = x[(i-0)*dof+j]-slope[(i-0)*dof+j]*hx/2;
213: }
214: (*ctx->physics.riemann)(ctx->physics.user,dof,uL,uR,ctx->flux,&maxspeed,ctx->xmin+hx*i,ctx->xmin,ctx->xmax);
215: cfl_idt = PetscMax(cfl_idt,PetscAbsScalar(maxspeed/hx)); /* Max allowable value of 1/Delta t */
216: if (i > xs) {
217: for (j=0; j<dof; j++) f[(i-1)*dof+j] -= ctx->flux[j]/hx;
218: }
219: if (i < xs+xm) {
220: for (j=0; j<dof; j++) f[i*dof+j] += ctx->flux[j]/hx;
221: }
222: }
223: if (i == len_slow) { /* slow parts can be changed based on a */
224: uL = &ctx->uLR[0];
225: uR = &ctx->uLR[dof];
226: for (j=0; j<dof; j++) {
227: uL[j] = x[(i-1)*dof+j]+slope[(i-1)*dof+j]*hx/2;
228: uR[j] = x[(i-0)*dof+j]-slope[(i-0)*dof+j]*hx/2;
229: }
230: (*ctx->physics.riemann)(ctx->physics.user,dof,uL,uR,ctx->flux,&maxspeed,ctx->xmin+hx*i,ctx->xmin,ctx->xmax);
231: cfl_idt = PetscMax(cfl_idt,PetscAbsScalar(maxspeed/hx)); /* Max allowable value of 1/Delta t */
232: if (i > xs) {
233: for (j=0; j<dof; j++) f[(i-1)*dof+j] -= ctx->flux[j]/hx;
234: }
235: }
236: }
237: DMDAVecRestoreArray(da,Xloc,&x);
238: VecRestoreArray(F,&f);
239: DMDARestoreArray(da,PETSC_TRUE,&slope);
240: DMRestoreLocalVector(da,&Xloc);
241: return(0);
242: }
244: /* --------------------------------- Finite Volume Solver for fast parts ----------------------------------- */
246: PetscErrorCode FVRHSFunctionfast(TS ts,PetscReal time,Vec X,Vec F,void *vctx)
247: {
248: FVCtx *ctx = (FVCtx*)vctx;
250: PetscInt i,j,k,Mx,dof,xs,xm,len_slow;
251: PetscReal hx,cfl_idt = 0;
252: PetscScalar *x,*f,*slope;
253: Vec Xloc;
254: DM da;
257: TSGetDM(ts,&da);
258: DMGetLocalVector(da,&Xloc);
259: DMDAGetInfo(da,0, &Mx,0,0, 0,0,0, &dof,0,0,0,0,0);
260: hx = (ctx->xmax-ctx->xmin)/Mx;
261: DMGlobalToLocalBegin(da,X,INSERT_VALUES,Xloc);
262: DMGlobalToLocalEnd(da,X,INSERT_VALUES,Xloc);
263: VecZeroEntries(F);
264: DMDAVecGetArray(da,Xloc,&x);
265: VecGetArray(F,&f);
266: DMDAGetArray(da,PETSC_TRUE,&slope);
267: DMDAGetCorners(da,&xs,0,0,&xm,0,0);
268: ISGetSize(ctx->iss,&len_slow);
270: if (ctx->bctype == FVBC_OUTFLOW) {
271: for (i=xs-2; i<0; i++) {
272: for (j=0; j<dof; j++) x[i*dof+j] = x[j];
273: }
274: for (i=Mx; i<xs+xm+2; i++) {
275: for (j=0; j<dof; j++) x[i*dof+j] = x[(xs+xm-1)*dof+j];
276: }
277: }
278: for (i=xs-1; i<xs+xm+1; i++) {
279: struct _LimitInfo info;
280: PetscScalar *cjmpL,*cjmpR;
281: if (i > len_slow-2) {
282: (*ctx->physics.characteristic)(ctx->physics.user,dof,&x[i*dof],ctx->R,ctx->Rinv,ctx->speeds,ctx->xmin+hx*i);
283: PetscArrayzero(ctx->cjmpLR,2*dof);
284: cjmpL = &ctx->cjmpLR[0];
285: cjmpR = &ctx->cjmpLR[dof];
286: for (j=0; j<dof; j++) {
287: PetscScalar jmpL,jmpR;
288: jmpL = x[(i+0)*dof+j]-x[(i-1)*dof+j];
289: jmpR = x[(i+1)*dof+j]-x[(i+0)*dof+j];
290: for (k=0; k<dof; k++) {
291: cjmpL[k] += ctx->Rinv[k+j*dof]*jmpL;
292: cjmpR[k] += ctx->Rinv[k+j*dof]*jmpR;
293: }
294: }
295: /* Apply limiter to the left and right characteristic jumps */
296: info.m = dof;
297: info.hx = hx;
298: (*ctx->limit)(&info,cjmpL,cjmpR,ctx->cslope);
299: for (j=0; j<dof; j++) ctx->cslope[j] /= hx; /* rescale to a slope */
300: for (j=0; j<dof; j++) {
301: PetscScalar tmp = 0;
302: for (k=0; k<dof; k++) tmp += ctx->R[j+k*dof]*ctx->cslope[k];
303: slope[i*dof+j] = tmp;
304: }
305: }
306: }
308: for (i=xs; i<xs+xm+1; i++) {
309: PetscReal maxspeed;
310: PetscScalar *uL,*uR;
311: if (i > len_slow) {
312: uL = &ctx->uLR[0];
313: uR = &ctx->uLR[dof];
314: for (j=0; j<dof; j++) {
315: uL[j] = x[(i-1)*dof+j]+slope[(i-1)*dof+j]*hx/2;
316: uR[j] = x[(i-0)*dof+j]-slope[(i-0)*dof+j]*hx/2;
317: }
318: (*ctx->physics.riemann)(ctx->physics.user,dof,uL,uR,ctx->flux,&maxspeed,ctx->xmin+hx*i,ctx->xmin,ctx->xmax);
319: cfl_idt = PetscMax(cfl_idt,PetscAbsScalar(maxspeed/hx)); /* Max allowable value of 1/Delta t */
320: if (i > xs) {
321: for (j=0; j<dof; j++) f[(i-len_slow-1)*dof+j] -= ctx->flux[j]/hx;
322: }
323: if (i < xs+xm) {
324: for (j=0; j<dof; j++) f[(i-len_slow)*dof+j] += ctx->flux[j]/hx;
325: }
326: }
327: if (i == len_slow) {
328: uL = &ctx->uLR[0];
329: uR = &ctx->uLR[dof];
330: for (j=0; j<dof; j++) {
331: uL[j] = x[(i-1)*dof+j]+slope[(i-1)*dof+j]*hx/2;
332: uR[j] = x[(i-0)*dof+j]-slope[(i-0)*dof+j]*hx/2;
333: }
334: (*ctx->physics.riemann)(ctx->physics.user,dof,uL,uR,ctx->flux,&maxspeed,ctx->xmin+hx*i,ctx->xmin,ctx->xmax);
335: cfl_idt = PetscMax(cfl_idt,PetscAbsScalar(maxspeed/hx)); /* Max allowable value of 1/Delta t */
336: if (i < xs+xm) {
337: for (j=0; j<dof; j++) f[(i-len_slow)*dof+j] += ctx->flux[j]/hx;
338: }
339: }
340: }
341: DMDAVecRestoreArray(da,Xloc,&x);
342: VecRestoreArray(F,&f);
343: DMDARestoreArray(da,PETSC_TRUE,&slope);
344: DMRestoreLocalVector(da,&Xloc);
345: return(0);
346: }
348: PetscErrorCode FVRHSFunctionslow2(TS ts,PetscReal time,Vec X,Vec F,void *vctx)
349: {
350: FVCtx *ctx = (FVCtx*)vctx;
352: PetscInt i,j,k,Mx,dof,xs,xm,len_slow1,len_slow2;
353: PetscReal hx,cfl_idt = 0;
354: PetscScalar *x,*f,*slope;
355: Vec Xloc;
356: DM da;
359: TSGetDM(ts,&da);
360: DMGetLocalVector(da,&Xloc);
361: DMDAGetInfo(da,0, &Mx,0,0, 0,0,0, &dof,0,0,0,0,0);
362: hx = (ctx->xmax-ctx->xmin)/Mx;
363: DMGlobalToLocalBegin(da,X,INSERT_VALUES,Xloc);
364: DMGlobalToLocalEnd(da,X,INSERT_VALUES,Xloc);
365: VecZeroEntries(F);
366: DMDAVecGetArray(da,Xloc,&x);
367: VecGetArray(F,&f);
368: DMDAGetArray(da,PETSC_TRUE,&slope);
369: DMDAGetCorners(da,&xs,0,0,&xm,0,0);
370: ISGetSize(ctx->iss,&len_slow1);
371: ISGetSize(ctx->iss2,&len_slow2);
372: if (ctx->bctype == FVBC_OUTFLOW) {
373: for (i=xs-2; i<0; i++) {
374: for (j=0; j<dof; j++) x[i*dof+j] = x[j];
375: }
376: for (i=Mx; i<xs+xm+2; i++) {
377: for (j=0; j<dof; j++) x[i*dof+j] = x[(xs+xm-1)*dof+j];
378: }
379: }
380: for (i=xs-1; i<xs+xm+1; i++) {
381: struct _LimitInfo info;
382: PetscScalar *cjmpL,*cjmpR;
383: if (i < len_slow1+len_slow2+1 && i > len_slow1-2) {
384: (*ctx->physics.characteristic)(ctx->physics.user,dof,&x[i*dof],ctx->R,ctx->Rinv,ctx->speeds,ctx->xmin+hx*i);
385: PetscArrayzero(ctx->cjmpLR,2*dof);
386: cjmpL = &ctx->cjmpLR[0];
387: cjmpR = &ctx->cjmpLR[dof];
388: for (j=0; j<dof; j++) {
389: PetscScalar jmpL,jmpR;
390: jmpL = x[(i+0)*dof+j]-x[(i-1)*dof+j];
391: jmpR = x[(i+1)*dof+j]-x[(i+0)*dof+j];
392: for (k=0; k<dof; k++) {
393: cjmpL[k] += ctx->Rinv[k+j*dof]*jmpL;
394: cjmpR[k] += ctx->Rinv[k+j*dof]*jmpR;
395: }
396: }
397: /* Apply limiter to the left and right characteristic jumps */
398: info.m = dof;
399: info.hx = hx;
400: (*ctx->limit)(&info,cjmpL,cjmpR,ctx->cslope);
401: for (j=0; j<dof; j++) ctx->cslope[j] /= hx; /* rescale to a slope */
402: for (j=0; j<dof; j++) {
403: PetscScalar tmp = 0;
404: for (k=0; k<dof; k++) tmp += ctx->R[j+k*dof]*ctx->cslope[k];
405: slope[i*dof+j] = tmp;
406: }
407: }
408: }
410: for (i=xs; i<xs+xm+1; i++) {
411: PetscReal maxspeed;
412: PetscScalar *uL,*uR;
413: if (i < len_slow1+len_slow2 && i > len_slow1) { /* slow parts can be changed based on a */
414: uL = &ctx->uLR[0];
415: uR = &ctx->uLR[dof];
416: for (j=0; j<dof; j++) {
417: uL[j] = x[(i-1)*dof+j]+slope[(i-1)*dof+j]*hx/2;
418: uR[j] = x[(i-0)*dof+j]-slope[(i-0)*dof+j]*hx/2;
419: }
420: (*ctx->physics.riemann)(ctx->physics.user,dof,uL,uR,ctx->flux,&maxspeed,ctx->xmin+hx*i,ctx->xmin,ctx->xmax);
421: cfl_idt = PetscMax(cfl_idt,PetscAbsScalar(maxspeed/hx)); /* Max allowable value of 1/Delta t */
422: if (i > xs) {
423: for (j=0; j<dof; j++) f[(i-len_slow1-1)*dof+j] -= ctx->flux[j]/hx;
424: }
425: if (i < xs+xm) {
426: for (j=0; j<dof; j++) f[(i-len_slow1)*dof+j] += ctx->flux[j]/hx;
427: }
428: }
429: if (i == len_slow1+len_slow2) { /* slow parts can be changed based on a */
430: uL = &ctx->uLR[0];
431: uR = &ctx->uLR[dof];
432: for (j=0; j<dof; j++) {
433: uL[j] = x[(i-1)*dof+j]+slope[(i-1)*dof+j]*hx/2;
434: uR[j] = x[(i-0)*dof+j]-slope[(i-0)*dof+j]*hx/2;
435: }
436: (*ctx->physics.riemann)(ctx->physics.user,dof,uL,uR,ctx->flux,&maxspeed,ctx->xmin+hx*i,ctx->xmin,ctx->xmax);
437: cfl_idt = PetscMax(cfl_idt,PetscAbsScalar(maxspeed/hx)); /* Max allowable value of 1/Delta t */
438: if (i > xs) {
439: for (j=0; j<dof; j++) f[(i-len_slow1-1)*dof+j] -= ctx->flux[j]/hx;
440: }
441: }
442: if (i == len_slow1) { /* slow parts can be changed based on a */
443: uL = &ctx->uLR[0];
444: uR = &ctx->uLR[dof];
445: for (j=0; j<dof; j++) {
446: uL[j] = x[(i-1)*dof+j]+slope[(i-1)*dof+j]*hx/2;
447: uR[j] = x[(i-0)*dof+j]-slope[(i-0)*dof+j]*hx/2;
448: }
449: (*ctx->physics.riemann)(ctx->physics.user,dof,uL,uR,ctx->flux,&maxspeed,ctx->xmin+hx*i,ctx->xmin,ctx->xmax);
450: cfl_idt = PetscMax(cfl_idt,PetscAbsScalar(maxspeed/hx)); /* Max allowable value of 1/Delta t */
451: if (i < xs+xm) {
452: for (j=0; j<dof; j++) f[(i-len_slow1)*dof+j] += ctx->flux[j]/hx;
453: }
454: }
455: }
457: DMDAVecRestoreArray(da,Xloc,&x);
458: VecRestoreArray(F,&f);
459: DMDARestoreArray(da,PETSC_TRUE,&slope);
460: DMRestoreLocalVector(da,&Xloc);
461: return(0);
462: }
464: PetscErrorCode FVRHSFunctionfast2(TS ts,PetscReal time,Vec X,Vec F,void *vctx)
465: {
466: FVCtx *ctx = (FVCtx*)vctx;
468: PetscInt i,j,k,Mx,dof,xs,xm,len_slow1,len_slow2;
469: PetscReal hx,cfl_idt = 0;
470: PetscScalar *x,*f,*slope;
471: Vec Xloc;
472: DM da;
475: TSGetDM(ts,&da);
476: DMGetLocalVector(da,&Xloc);
477: DMDAGetInfo(da,0, &Mx,0,0, 0,0,0, &dof,0,0,0,0,0);
478: hx = (ctx->xmax-ctx->xmin)/Mx;
479: DMGlobalToLocalBegin(da,X,INSERT_VALUES,Xloc);
480: DMGlobalToLocalEnd(da,X,INSERT_VALUES,Xloc);
481: VecZeroEntries(F);
482: DMDAVecGetArray(da,Xloc,&x);
483: VecGetArray(F,&f);
484: DMDAGetArray(da,PETSC_TRUE,&slope);
485: DMDAGetCorners(da,&xs,0,0,&xm,0,0);
486: ISGetSize(ctx->iss,&len_slow1);
487: ISGetSize(ctx->iss2,&len_slow2);
489: if (ctx->bctype == FVBC_OUTFLOW) {
490: for (i=xs-2; i<0; i++) {
491: for (j=0; j<dof; j++) x[i*dof+j] = x[j];
492: }
493: for (i=Mx; i<xs+xm+2; i++) {
494: for (j=0; j<dof; j++) x[i*dof+j] = x[(xs+xm-1)*dof+j];
495: }
496: }
497: for (i=xs-1; i<xs+xm+1; i++) {
498: struct _LimitInfo info;
499: PetscScalar *cjmpL,*cjmpR;
500: if (i > len_slow1+len_slow2-2) {
501: (*ctx->physics.characteristic)(ctx->physics.user,dof,&x[i*dof],ctx->R,ctx->Rinv,ctx->speeds,ctx->xmin+hx*i);
502: PetscArrayzero(ctx->cjmpLR,2*dof);
503: cjmpL = &ctx->cjmpLR[0];
504: cjmpR = &ctx->cjmpLR[dof];
505: for (j=0; j<dof; j++) {
506: PetscScalar jmpL,jmpR;
507: jmpL = x[(i+0)*dof+j]-x[(i-1)*dof+j];
508: jmpR = x[(i+1)*dof+j]-x[(i+0)*dof+j];
509: for (k=0; k<dof; k++) {
510: cjmpL[k] += ctx->Rinv[k+j*dof]*jmpL;
511: cjmpR[k] += ctx->Rinv[k+j*dof]*jmpR;
512: }
513: }
514: /* Apply limiter to the left and right characteristic jumps */
515: info.m = dof;
516: info.hx = hx;
517: (*ctx->limit)(&info,cjmpL,cjmpR,ctx->cslope);
518: for (j=0; j<dof; j++) ctx->cslope[j] /= hx; /* rescale to a slope */
519: for (j=0; j<dof; j++) {
520: PetscScalar tmp = 0;
521: for (k=0; k<dof; k++) tmp += ctx->R[j+k*dof]*ctx->cslope[k];
522: slope[i*dof+j] = tmp;
523: }
524: }
525: }
527: for (i=xs; i<xs+xm+1; i++) {
528: PetscReal maxspeed;
529: PetscScalar *uL,*uR;
530: if (i > len_slow1+len_slow2) {
531: uL = &ctx->uLR[0];
532: uR = &ctx->uLR[dof];
533: for (j=0; j<dof; j++) {
534: uL[j] = x[(i-1)*dof+j]+slope[(i-1)*dof+j]*hx/2;
535: uR[j] = x[(i-0)*dof+j]-slope[(i-0)*dof+j]*hx/2;
536: }
537: (*ctx->physics.riemann)(ctx->physics.user,dof,uL,uR,ctx->flux,&maxspeed,ctx->xmin+hx*i,ctx->xmin,ctx->xmax);
538: cfl_idt = PetscMax(cfl_idt,PetscAbsScalar(maxspeed/hx)); /* Max allowable value of 1/Delta t */
539: if (i > xs) {
540: for (j=0; j<dof; j++) f[(i-len_slow1-len_slow2-1)*dof+j] -= ctx->flux[j]/hx;
541: }
542: if (i < xs+xm) {
543: for (j=0; j<dof; j++) f[(i-len_slow1-len_slow2)*dof+j] += ctx->flux[j]/hx;
544: }
545: }
546: if (i == len_slow1+len_slow2) {
547: uL = &ctx->uLR[0];
548: uR = &ctx->uLR[dof];
549: for (j=0; j<dof; j++) {
550: uL[j] = x[(i-1)*dof+j]+slope[(i-1)*dof+j]*hx/2;
551: uR[j] = x[(i-0)*dof+j]-slope[(i-0)*dof+j]*hx/2;
552: }
553: (*ctx->physics.riemann)(ctx->physics.user,dof,uL,uR,ctx->flux,&maxspeed,ctx->xmin+hx*i,ctx->xmin,ctx->xmax);
554: cfl_idt = PetscMax(cfl_idt,PetscAbsScalar(maxspeed/hx)); /* Max allowable value of 1/Delta t */
555: if (i < xs+xm) {
556: for (j=0; j<dof; j++) f[(i-len_slow1-len_slow2)*dof+j] += ctx->flux[j]/hx;
557: }
558: }
559: }
560: DMDAVecRestoreArray(da,Xloc,&x);
561: VecRestoreArray(F,&f);
562: DMDARestoreArray(da,PETSC_TRUE,&slope);
563: DMRestoreLocalVector(da,&Xloc);
564: return(0);
565: }
567: int main(int argc,char *argv[])
568: {
569: char lname[256] = "mc",physname[256] = "advect",final_fname[256] = "solution.m";
570: PetscFunctionList limiters = 0,physics = 0;
571: MPI_Comm comm;
572: TS ts;
573: DM da;
574: Vec X,X0,R;
575: FVCtx ctx;
576: PetscInt i,k,dof,xs,xm,Mx,draw = 0,*index_slow,*index_fast,islow = 0,ifast = 0;
577: PetscBool view_final = PETSC_FALSE;
578: PetscReal ptime;
579: PetscErrorCode ierr;
581: PetscInitialize(&argc,&argv,0,help);if (ierr) return ierr;
582: comm = PETSC_COMM_WORLD;
583: PetscMemzero(&ctx,sizeof(ctx));
585: /* Register limiters to be available on the command line */
586: PetscFunctionListAdd(&limiters,"upwind" ,Limit_Upwind);
587: PetscFunctionListAdd(&limiters,"lax-wendroff" ,Limit_LaxWendroff);
588: PetscFunctionListAdd(&limiters,"beam-warming" ,Limit_BeamWarming);
589: PetscFunctionListAdd(&limiters,"fromm" ,Limit_Fromm);
590: PetscFunctionListAdd(&limiters,"minmod" ,Limit_Minmod);
591: PetscFunctionListAdd(&limiters,"superbee" ,Limit_Superbee);
592: PetscFunctionListAdd(&limiters,"mc" ,Limit_MC);
593: PetscFunctionListAdd(&limiters,"vanleer" ,Limit_VanLeer);
594: PetscFunctionListAdd(&limiters,"vanalbada" ,Limit_VanAlbada);
595: PetscFunctionListAdd(&limiters,"vanalbadatvd" ,Limit_VanAlbadaTVD);
596: PetscFunctionListAdd(&limiters,"koren" ,Limit_Koren);
597: PetscFunctionListAdd(&limiters,"korensym" ,Limit_KorenSym);
598: PetscFunctionListAdd(&limiters,"koren3" ,Limit_Koren3);
599: PetscFunctionListAdd(&limiters,"cada-torrilhon2" ,Limit_CadaTorrilhon2);
600: PetscFunctionListAdd(&limiters,"cada-torrilhon3-r0p1",Limit_CadaTorrilhon3R0p1);
601: PetscFunctionListAdd(&limiters,"cada-torrilhon3-r1" ,Limit_CadaTorrilhon3R1);
602: PetscFunctionListAdd(&limiters,"cada-torrilhon3-r10" ,Limit_CadaTorrilhon3R10);
603: PetscFunctionListAdd(&limiters,"cada-torrilhon3-r100",Limit_CadaTorrilhon3R100);
605: /* Register physical models to be available on the command line */
606: PetscFunctionListAdd(&physics,"advect" ,PhysicsCreate_Advect);
608: ctx.comm = comm;
609: ctx.cfl = 0.9;
610: ctx.bctype = FVBC_PERIODIC;
611: ctx.xmin = -1.0;
612: ctx.xmax = 1.0;
613: PetscOptionsBegin(comm,NULL,"Finite Volume solver options","");
614: PetscOptionsReal("-xmin","X min","",ctx.xmin,&ctx.xmin,NULL);
615: PetscOptionsReal("-xmax","X max","",ctx.xmax,&ctx.xmax,NULL);
616: PetscOptionsFList("-limit","Name of flux limiter to use","",limiters,lname,lname,sizeof(lname),NULL);
617: PetscOptionsInt("-draw","Draw solution vector, bitwise OR of (1=initial,2=final,4=final error)","",draw,&draw,NULL);
618: PetscOptionsString("-view_final","Write final solution in ASCII MATLAB format to given file name","",final_fname,final_fname,sizeof(final_fname),&view_final);
619: PetscOptionsInt("-initial","Initial condition (depends on the physics)","",ctx.initial,&ctx.initial,NULL);
620: PetscOptionsBool("-simulation","Compare errors with reference solution","",ctx.simulation,&ctx.simulation,NULL);
621: PetscOptionsReal("-cfl","CFL number to time step at","",ctx.cfl,&ctx.cfl,NULL);
622: PetscOptionsEnum("-bc_type","Boundary condition","",FVBCTypes,(PetscEnum)ctx.bctype,(PetscEnum*)&ctx.bctype,NULL);
623: PetscOptionsBool("-recursive_split","Split the domain recursively","",ctx.recursive,&ctx.recursive,NULL);
624: PetscOptionsEnd();
626: /* Choose the limiter from the list of registered limiters */
627: PetscFunctionListFind(limiters,lname,&ctx.limit);
628: if (!ctx.limit) SETERRQ1(PETSC_COMM_SELF,1,"Limiter '%s' not found",lname);
630: /* Choose the physics from the list of registered models */
631: {
632: PetscErrorCode (*r)(FVCtx*);
633: PetscFunctionListFind(physics,physname,&r);
634: if (!r) SETERRQ1(PETSC_COMM_SELF,1,"Physics '%s' not found",physname);
635: /* Create the physics, will set the number of fields and their names */
636: (*r)(&ctx);
637: }
639: /* Create a DMDA to manage the parallel grid */
640: DMDACreate1d(comm,DM_BOUNDARY_PERIODIC,50,ctx.physics.dof,2,NULL,&da);
641: DMSetFromOptions(da);
642: DMSetUp(da);
643: /* Inform the DMDA of the field names provided by the physics. */
644: /* The names will be shown in the title bars when run with -ts_monitor_draw_solution */
645: for (i=0; i<ctx.physics.dof; i++) {
646: DMDASetFieldName(da,i,ctx.physics.fieldname[i]);
647: }
648: DMDAGetInfo(da,0, &Mx,0,0, 0,0,0, &dof,0,0,0,0,0);
649: DMDAGetCorners(da,&xs,0,0,&xm,0,0);
651: /* Set coordinates of cell centers */
652: DMDASetUniformCoordinates(da,ctx.xmin+0.5*(ctx.xmax-ctx.xmin)/Mx,ctx.xmax+0.5*(ctx.xmax-ctx.xmin)/Mx,0,0,0,0);
654: /* Allocate work space for the Finite Volume solver (so it doesn't have to be reallocated on each function evaluation) */
655: PetscMalloc4(dof*dof,&ctx.R,dof*dof,&ctx.Rinv,2*dof,&ctx.cjmpLR,1*dof,&ctx.cslope);
656: PetscMalloc3(2*dof,&ctx.uLR,dof,&ctx.flux,dof,&ctx.speeds);
658: /* Create a vector to store the solution and to save the initial state */
659: DMCreateGlobalVector(da,&X);
660: VecDuplicate(X,&X0);
661: VecDuplicate(X,&R);
663: /*-------------------------------- create index for slow parts and fast parts ----------------------------------------*/
664: PetscMalloc1(xm*dof,&index_slow);
665: PetscMalloc1(xm*dof,&index_fast);
666: for (i=xs; i<xs+xm; i++) {
667: if (ctx.xmin+i*(ctx.xmax-ctx.xmin)/(PetscReal)Mx+0.5*(ctx.xmax-ctx.xmin)/(PetscReal)Mx < 0)
668: for (k=0; k<dof; k++) index_slow[islow++] = i*dof+k;
669: else
670: for (k=0; k<dof; k++) index_fast[ifast++] = i*dof+k;
671: } /* this step need to be changed based on discontinuous point of a */
672: ISCreateGeneral(PETSC_COMM_WORLD,islow,index_slow,PETSC_COPY_VALUES,&ctx.iss);
673: ISCreateGeneral(PETSC_COMM_WORLD,ifast,index_fast,PETSC_COPY_VALUES,&ctx.isf);
675: /* Create a time-stepping object */
676: TSCreate(comm,&ts);
677: TSSetDM(ts,da);
678: TSSetRHSFunction(ts,R,FVRHSFunction,&ctx);
679: TSRHSSplitSetIS(ts,"slow",ctx.iss);
680: TSRHSSplitSetIS(ts,"fast",ctx.isf);
681: TSRHSSplitSetRHSFunction(ts,"slow",NULL,FVRHSFunctionslow,&ctx);
682: TSRHSSplitSetRHSFunction(ts,"fast",NULL,FVRHSFunctionfast,&ctx);
684: if (ctx.recursive) {
685: TS subts;
686: islow = 0;
687: ifast = 0;
688: for (i=xs; i<xs+xm; i++) {
689: PetscReal coord = ctx.xmin+i*(ctx.xmax-ctx.xmin)/(PetscReal)Mx+0.5*(ctx.xmax-ctx.xmin)/(PetscReal)Mx;
690: if (coord >= 0 && coord < ctx.xmin+(ctx.xmax-ctx.xmin)*3/4.)
691: for (k=0; k<dof; k++) index_slow[islow++] = i*dof+k;
692: if (coord >= ctx.xmin+(ctx.xmax-ctx.xmin)*3/4.)
693: for (k=0; k<dof; k++) index_fast[ifast++] = i*dof+k;
694: } /* this step need to be changed based on discontinuous point of a */
695: ISCreateGeneral(PETSC_COMM_WORLD,islow,index_slow,PETSC_COPY_VALUES,&ctx.iss2);
696: ISCreateGeneral(PETSC_COMM_WORLD,ifast,index_fast,PETSC_COPY_VALUES,&ctx.isf2);
698: TSRHSSplitGetSubTS(ts,"fast",&subts);
699: TSRHSSplitSetIS(subts,"slow",ctx.iss2);
700: TSRHSSplitSetIS(subts,"fast",ctx.isf2);
701: TSRHSSplitSetRHSFunction(subts,"slow",NULL,FVRHSFunctionslow2,&ctx);
702: TSRHSSplitSetRHSFunction(subts,"fast",NULL,FVRHSFunctionfast2,&ctx);
703: }
705: /*TSSetType(ts,TSSSP);*/
706: TSSetType(ts,TSMPRK);
707: TSSetMaxTime(ts,10);
708: TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);
710: /* Compute initial conditions and starting time step */
711: FVSample(&ctx,da,0,X0);
712: FVRHSFunction(ts,0,X0,X,(void*)&ctx); /* Initial function evaluation, only used to determine max speed */
713: VecCopy(X0,X); /* The function value was not used so we set X=X0 again */
714: TSSetTimeStep(ts,ctx.cfl/ctx.cfl_idt);
715: TSSetFromOptions(ts); /* Take runtime options */
716: SolutionStatsView(da,X,PETSC_VIEWER_STDOUT_WORLD);
717: {
718: PetscInt steps;
719: PetscScalar mass_initial,mass_final,mass_difference;
721: TSSolve(ts,X);
722: TSGetSolveTime(ts,&ptime);
723: TSGetStepNumber(ts,&steps);
724: PetscPrintf(comm,"Final time %g, steps %D\n",(double)ptime,steps);
725: /* calculate the total mass at initial time and final time */
726: mass_initial = 0.0;
727: mass_final = 0.0;
728: VecSum(X0,&mass_initial);
729: VecSum(X,&mass_final);
730: mass_difference = (ctx.xmax-ctx.xmin)/(PetscScalar)Mx*(mass_final - mass_initial);
731: PetscPrintf(comm,"Mass difference %g\n",(double)mass_difference);
732: if (ctx.simulation) {
733: PetscViewer fd;
734: char filename[PETSC_MAX_PATH_LEN] = "binaryoutput";
735: Vec XR;
736: PetscReal nrm1,nrmsup;
737: PetscBool flg;
738: PetscOptionsGetString(NULL,NULL,"-f",filename,sizeof(filename),&flg);
739: if (!flg) SETERRQ(PETSC_COMM_WORLD,1,"Must indicate binary file with the -f option");
740: PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_READ,&fd);
741: VecDuplicate(X0,&XR);
742: VecLoad(XR,fd);
743: PetscViewerDestroy(&fd);
744: VecAYPX(XR,-1,X);
745: VecNorm(XR,NORM_1,&nrm1);
746: VecNorm(XR,NORM_INFINITY,&nrmsup);
747: nrm1 /= Mx;
748: PetscPrintf(comm,"Error ||x-x_e||_1 %g ||x-x_e||_sup %g\n",(double)nrm1,(double)nrmsup);
749: VecDestroy(&XR);
750: }
751: }
753: SolutionStatsView(da,X,PETSC_VIEWER_STDOUT_WORLD);
754: if (draw & 0x1) {VecView(X0,PETSC_VIEWER_DRAW_WORLD);}
755: if (draw & 0x2) {VecView(X,PETSC_VIEWER_DRAW_WORLD);}
756: if (draw & 0x4) {
757: Vec Y;
758: VecDuplicate(X,&Y);
759: FVSample(&ctx,da,ptime,Y);
760: VecAYPX(Y,-1,X);
761: VecView(Y,PETSC_VIEWER_DRAW_WORLD);
762: VecDestroy(&Y);
763: }
765: if (view_final) {
766: PetscViewer viewer;
767: PetscViewerASCIIOpen(PETSC_COMM_WORLD,final_fname,&viewer);
768: PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_MATLAB);
769: VecView(X,viewer);
770: PetscViewerPopFormat(viewer);
771: PetscViewerDestroy(&viewer);
772: }
774: /* Clean up */
775: (*ctx.physics.destroy)(ctx.physics.user);
776: for (i=0; i<ctx.physics.dof; i++) {PetscFree(ctx.physics.fieldname[i]);}
777: PetscFree4(ctx.R,ctx.Rinv,ctx.cjmpLR,ctx.cslope);
778: PetscFree3(ctx.uLR,ctx.flux,ctx.speeds);
779: ISDestroy(&ctx.iss);
780: ISDestroy(&ctx.isf);
781: ISDestroy(&ctx.iss2);
782: ISDestroy(&ctx.isf2);
783: VecDestroy(&X);
784: VecDestroy(&X0);
785: VecDestroy(&R);
786: DMDestroy(&da);
787: TSDestroy(&ts);
788: PetscFree(index_slow);
789: PetscFree(index_fast);
790: PetscFunctionListDestroy(&limiters);
791: PetscFunctionListDestroy(&physics);
792: PetscFinalize();
793: return ierr;
794: }
796: /*TEST
797: build:
798: requires: !complex c99
799: depends: finitevolume1d.c
801: test:
802: 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
804: test:
805: suffix: 2
806: 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
807: output_file: output/ex5_1.out
809: test:
810: suffix: 3
811: 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
813: test:
814: suffix: 4
815: 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
816: output_file: output/ex5_3.out
818: test:
819: suffix: 5
820: 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
822: test:
823: suffix: 6
824: 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
825: TEST*/