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
5
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 */
12
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";
23
24 #include <petscts.h>
25 #include <petscdm.h>
26 #include <petscdmda.h>
27 #include <petscdraw.h>
28 #include <petsc/private/tsimpl.h>
29
30 #include <petsc/private/kernels/blockinvert.h> /* For the Kernel_*_gets_* stuff for BAIJ */
31
32 #include "finitevolume1d.h"
33
RangeMod(PetscReal a,PetscReal xmin,PetscReal xmax)34 PETSC_STATIC_INLINE PetscReal RangeMod(PetscReal a,PetscReal xmin,PetscReal xmax) { PetscReal range = xmax-xmin; return xmin +PetscFmodReal(range+PetscFmodReal(a,range),range); }
35
36 /* --------------------------------- Advection ----------------------------------- */
37
38 typedef struct {
39 PetscReal a[2]; /* advective velocity */
40 } AdvectCtx;
41
PhysicsRiemann_Advect(void * vctx,PetscInt m,const PetscScalar * uL,const PetscScalar * uR,PetscScalar * flux,PetscReal * maxspeed,PetscReal x,PetscReal xmin,PetscReal xmax)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;
46
47 PetscFunctionBeginUser;
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 PetscFunctionReturn(0);
54 }
55
PhysicsCharacteristic_Advect(void * vctx,PetscInt m,const PetscScalar * u,PetscScalar * X,PetscScalar * Xi,PetscReal * speeds,PetscReal x)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;
59
60 PetscFunctionBeginUser;
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 PetscFunctionReturn(0);
66 }
67
PhysicsSample_Advect(void * vctx,PetscInt initial,FVBCType bctype,PetscReal xmin,PetscReal xmax,PetscReal t,PetscReal x,PetscReal * u)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;
72
73 PetscFunctionBeginUser;
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 PetscFunctionReturn(0);
111 }
112
PhysicsCreate_Advect(FVCtx * ctx)113 static PetscErrorCode PhysicsCreate_Advect(FVCtx *ctx)
114 {
115 PetscErrorCode ierr;
116 AdvectCtx *user;
117
118 PetscFunctionBeginUser;
119 ierr = PetscNew(&user);CHKERRQ(ierr);
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 ierr = PetscStrallocpy("u",&ctx->physics.fieldname[0]);CHKERRQ(ierr);
127 user->a[0] = 1;
128 user->a[1] = 1;
129 ierr = PetscOptionsBegin(ctx->comm,ctx->prefix,"Options for advection","");CHKERRQ(ierr);
130 {
131 ierr = PetscOptionsReal("-physics_advect_a1","Speed1","",user->a[0],&user->a[0],NULL);CHKERRQ(ierr);
132 ierr = PetscOptionsReal("-physics_advect_a2","Speed2","",user->a[1],&user->a[1],NULL);CHKERRQ(ierr);
133 }
134 ierr = PetscOptionsEnd();CHKERRQ(ierr);
135 PetscFunctionReturn(0);
136 }
137
138 /* --------------------------------- Finite Volume Solver for slow parts ----------------------------------- */
139
FVRHSFunctionslow(TS ts,PetscReal time,Vec X,Vec F,void * vctx)140 PetscErrorCode FVRHSFunctionslow(TS ts,PetscReal time,Vec X,Vec F,void *vctx)
141 {
142 FVCtx *ctx = (FVCtx*)vctx;
143 PetscErrorCode ierr;
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;
149
150 PetscFunctionBeginUser;
151 ierr = TSGetDM(ts,&da);CHKERRQ(ierr);
152 ierr = DMGetLocalVector(da,&Xloc);CHKERRQ(ierr);
153 ierr = DMDAGetInfo(da,0, &Mx,0,0, 0,0,0, &dof,0,0,0,0,0);CHKERRQ(ierr);
154 hx = (ctx->xmax-ctx->xmin)/Mx;
155 ierr = DMGlobalToLocalBegin(da,X,INSERT_VALUES,Xloc);CHKERRQ(ierr);
156 ierr = DMGlobalToLocalEnd(da,X,INSERT_VALUES,Xloc);CHKERRQ(ierr);
157 ierr = VecZeroEntries(F);CHKERRQ(ierr);
158 ierr = DMDAVecGetArray(da,Xloc,&x);CHKERRQ(ierr);
159 ierr = VecGetArray(F,&f);CHKERRQ(ierr);
160 ierr = DMDAGetArray(da,PETSC_TRUE,&slope);CHKERRQ(ierr);
161 ierr = DMDAGetCorners(da,&xs,0,0,&xm,0,0);CHKERRQ(ierr);
162 ierr = ISGetSize(ctx->iss,&len_slow);CHKERRQ(ierr);
163
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 ierr = (*ctx->physics.characteristic)(ctx->physics.user,dof,&x[i*dof],ctx->R,ctx->Rinv,ctx->speeds,ctx->xmin+hx*i);CHKERRQ(ierr);
178 /* Evaluate jumps across interfaces (i-1, i) and (i, i+1), put in characteristic basis */
179 ierr = PetscArrayzero(ctx->cjmpLR,2*dof);CHKERRQ(ierr);
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 }
203
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 ierr = (*ctx->physics.riemann)(ctx->physics.user,dof,uL,uR,ctx->flux,&maxspeed,ctx->xmin+hx*i,ctx->xmin,ctx->xmax);CHKERRQ(ierr);
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 ierr = (*ctx->physics.riemann)(ctx->physics.user,dof,uL,uR,ctx->flux,&maxspeed,ctx->xmin+hx*i,ctx->xmin,ctx->xmax);CHKERRQ(ierr);
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 ierr = DMDAVecRestoreArray(da,Xloc,&x);CHKERRQ(ierr);
238 ierr = VecRestoreArray(F,&f);CHKERRQ(ierr);
239 ierr = DMDARestoreArray(da,PETSC_TRUE,&slope);CHKERRQ(ierr);
240 ierr = DMRestoreLocalVector(da,&Xloc);CHKERRQ(ierr);
241 PetscFunctionReturn(0);
242 }
243
244 /* --------------------------------- Finite Volume Solver for fast parts ----------------------------------- */
245
FVRHSFunctionfast(TS ts,PetscReal time,Vec X,Vec F,void * vctx)246 PetscErrorCode FVRHSFunctionfast(TS ts,PetscReal time,Vec X,Vec F,void *vctx)
247 {
248 FVCtx *ctx = (FVCtx*)vctx;
249 PetscErrorCode ierr;
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;
255
256 PetscFunctionBeginUser;
257 ierr = TSGetDM(ts,&da);CHKERRQ(ierr);
258 ierr = DMGetLocalVector(da,&Xloc);CHKERRQ(ierr);
259 ierr = DMDAGetInfo(da,0, &Mx,0,0, 0,0,0, &dof,0,0,0,0,0);CHKERRQ(ierr);
260 hx = (ctx->xmax-ctx->xmin)/Mx;
261 ierr = DMGlobalToLocalBegin(da,X,INSERT_VALUES,Xloc);CHKERRQ(ierr);
262 ierr = DMGlobalToLocalEnd(da,X,INSERT_VALUES,Xloc);CHKERRQ(ierr);
263 ierr = VecZeroEntries(F);CHKERRQ(ierr);
264 ierr = DMDAVecGetArray(da,Xloc,&x);CHKERRQ(ierr);
265 ierr = VecGetArray(F,&f);CHKERRQ(ierr);
266 ierr = DMDAGetArray(da,PETSC_TRUE,&slope);CHKERRQ(ierr);
267 ierr = DMDAGetCorners(da,&xs,0,0,&xm,0,0);CHKERRQ(ierr);
268 ierr = ISGetSize(ctx->iss,&len_slow);CHKERRQ(ierr);
269
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 ierr = (*ctx->physics.characteristic)(ctx->physics.user,dof,&x[i*dof],ctx->R,ctx->Rinv,ctx->speeds,ctx->xmin+hx*i);CHKERRQ(ierr);
283 ierr = PetscArrayzero(ctx->cjmpLR,2*dof);CHKERRQ(ierr);
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 }
307
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 ierr = (*ctx->physics.riemann)(ctx->physics.user,dof,uL,uR,ctx->flux,&maxspeed,ctx->xmin+hx*i,ctx->xmin,ctx->xmax);CHKERRQ(ierr);
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 ierr = (*ctx->physics.riemann)(ctx->physics.user,dof,uL,uR,ctx->flux,&maxspeed,ctx->xmin+hx*i,ctx->xmin,ctx->xmax);CHKERRQ(ierr);
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 ierr = DMDAVecRestoreArray(da,Xloc,&x);CHKERRQ(ierr);
342 ierr = VecRestoreArray(F,&f);CHKERRQ(ierr);
343 ierr = DMDARestoreArray(da,PETSC_TRUE,&slope);CHKERRQ(ierr);
344 ierr = DMRestoreLocalVector(da,&Xloc);CHKERRQ(ierr);
345 PetscFunctionReturn(0);
346 }
347
FVRHSFunctionslow2(TS ts,PetscReal time,Vec X,Vec F,void * vctx)348 PetscErrorCode FVRHSFunctionslow2(TS ts,PetscReal time,Vec X,Vec F,void *vctx)
349 {
350 FVCtx *ctx = (FVCtx*)vctx;
351 PetscErrorCode ierr;
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;
357
358 PetscFunctionBeginUser;
359 ierr = TSGetDM(ts,&da);CHKERRQ(ierr);
360 ierr = DMGetLocalVector(da,&Xloc);CHKERRQ(ierr);
361 ierr = DMDAGetInfo(da,0, &Mx,0,0, 0,0,0, &dof,0,0,0,0,0);CHKERRQ(ierr);
362 hx = (ctx->xmax-ctx->xmin)/Mx;
363 ierr = DMGlobalToLocalBegin(da,X,INSERT_VALUES,Xloc);CHKERRQ(ierr);
364 ierr = DMGlobalToLocalEnd(da,X,INSERT_VALUES,Xloc);CHKERRQ(ierr);
365 ierr = VecZeroEntries(F);CHKERRQ(ierr);
366 ierr = DMDAVecGetArray(da,Xloc,&x);CHKERRQ(ierr);
367 ierr = VecGetArray(F,&f);CHKERRQ(ierr);
368 ierr = DMDAGetArray(da,PETSC_TRUE,&slope);CHKERRQ(ierr);
369 ierr = DMDAGetCorners(da,&xs,0,0,&xm,0,0);CHKERRQ(ierr);
370 ierr = ISGetSize(ctx->iss,&len_slow1);CHKERRQ(ierr);
371 ierr = ISGetSize(ctx->iss2,&len_slow2);CHKERRQ(ierr);
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 ierr = (*ctx->physics.characteristic)(ctx->physics.user,dof,&x[i*dof],ctx->R,ctx->Rinv,ctx->speeds,ctx->xmin+hx*i);CHKERRQ(ierr);
385 ierr = PetscArrayzero(ctx->cjmpLR,2*dof);CHKERRQ(ierr);
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 }
409
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 ierr = (*ctx->physics.riemann)(ctx->physics.user,dof,uL,uR,ctx->flux,&maxspeed,ctx->xmin+hx*i,ctx->xmin,ctx->xmax);CHKERRQ(ierr);
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 ierr = (*ctx->physics.riemann)(ctx->physics.user,dof,uL,uR,ctx->flux,&maxspeed,ctx->xmin+hx*i,ctx->xmin,ctx->xmax);CHKERRQ(ierr);
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 ierr = (*ctx->physics.riemann)(ctx->physics.user,dof,uL,uR,ctx->flux,&maxspeed,ctx->xmin+hx*i,ctx->xmin,ctx->xmax);CHKERRQ(ierr);
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 }
456
457 ierr = DMDAVecRestoreArray(da,Xloc,&x);CHKERRQ(ierr);
458 ierr = VecRestoreArray(F,&f);CHKERRQ(ierr);
459 ierr = DMDARestoreArray(da,PETSC_TRUE,&slope);CHKERRQ(ierr);
460 ierr = DMRestoreLocalVector(da,&Xloc);CHKERRQ(ierr);
461 PetscFunctionReturn(0);
462 }
463
FVRHSFunctionfast2(TS ts,PetscReal time,Vec X,Vec F,void * vctx)464 PetscErrorCode FVRHSFunctionfast2(TS ts,PetscReal time,Vec X,Vec F,void *vctx)
465 {
466 FVCtx *ctx = (FVCtx*)vctx;
467 PetscErrorCode ierr;
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;
473
474 PetscFunctionBeginUser;
475 ierr = TSGetDM(ts,&da);CHKERRQ(ierr);
476 ierr = DMGetLocalVector(da,&Xloc);CHKERRQ(ierr);
477 ierr = DMDAGetInfo(da,0, &Mx,0,0, 0,0,0, &dof,0,0,0,0,0);CHKERRQ(ierr);
478 hx = (ctx->xmax-ctx->xmin)/Mx;
479 ierr = DMGlobalToLocalBegin(da,X,INSERT_VALUES,Xloc);CHKERRQ(ierr);
480 ierr = DMGlobalToLocalEnd(da,X,INSERT_VALUES,Xloc);CHKERRQ(ierr);
481 ierr = VecZeroEntries(F);CHKERRQ(ierr);
482 ierr = DMDAVecGetArray(da,Xloc,&x);CHKERRQ(ierr);
483 ierr = VecGetArray(F,&f);CHKERRQ(ierr);
484 ierr = DMDAGetArray(da,PETSC_TRUE,&slope);CHKERRQ(ierr);
485 ierr = DMDAGetCorners(da,&xs,0,0,&xm,0,0);CHKERRQ(ierr);
486 ierr = ISGetSize(ctx->iss,&len_slow1);CHKERRQ(ierr);
487 ierr = ISGetSize(ctx->iss2,&len_slow2);CHKERRQ(ierr);
488
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 ierr = (*ctx->physics.characteristic)(ctx->physics.user,dof,&x[i*dof],ctx->R,ctx->Rinv,ctx->speeds,ctx->xmin+hx*i);CHKERRQ(ierr);
502 ierr = PetscArrayzero(ctx->cjmpLR,2*dof);CHKERRQ(ierr);
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 }
526
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 ierr = (*ctx->physics.riemann)(ctx->physics.user,dof,uL,uR,ctx->flux,&maxspeed,ctx->xmin+hx*i,ctx->xmin,ctx->xmax);CHKERRQ(ierr);
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 ierr = (*ctx->physics.riemann)(ctx->physics.user,dof,uL,uR,ctx->flux,&maxspeed,ctx->xmin+hx*i,ctx->xmin,ctx->xmax);CHKERRQ(ierr);
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 ierr = DMDAVecRestoreArray(da,Xloc,&x);CHKERRQ(ierr);
561 ierr = VecRestoreArray(F,&f);CHKERRQ(ierr);
562 ierr = DMDARestoreArray(da,PETSC_TRUE,&slope);CHKERRQ(ierr);
563 ierr = DMRestoreLocalVector(da,&Xloc);CHKERRQ(ierr);
564 PetscFunctionReturn(0);
565 }
566
main(int argc,char * argv[])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;
580
581 ierr = PetscInitialize(&argc,&argv,0,help);if (ierr) return ierr;
582 comm = PETSC_COMM_WORLD;
583 ierr = PetscMemzero(&ctx,sizeof(ctx));CHKERRQ(ierr);
584
585 /* Register limiters to be available on the command line */
586 ierr = PetscFunctionListAdd(&limiters,"upwind" ,Limit_Upwind);CHKERRQ(ierr);
587 ierr = PetscFunctionListAdd(&limiters,"lax-wendroff" ,Limit_LaxWendroff);CHKERRQ(ierr);
588 ierr = PetscFunctionListAdd(&limiters,"beam-warming" ,Limit_BeamWarming);CHKERRQ(ierr);
589 ierr = PetscFunctionListAdd(&limiters,"fromm" ,Limit_Fromm);CHKERRQ(ierr);
590 ierr = PetscFunctionListAdd(&limiters,"minmod" ,Limit_Minmod);CHKERRQ(ierr);
591 ierr = PetscFunctionListAdd(&limiters,"superbee" ,Limit_Superbee);CHKERRQ(ierr);
592 ierr = PetscFunctionListAdd(&limiters,"mc" ,Limit_MC);CHKERRQ(ierr);
593 ierr = PetscFunctionListAdd(&limiters,"vanleer" ,Limit_VanLeer);CHKERRQ(ierr);
594 ierr = PetscFunctionListAdd(&limiters,"vanalbada" ,Limit_VanAlbada);CHKERRQ(ierr);
595 ierr = PetscFunctionListAdd(&limiters,"vanalbadatvd" ,Limit_VanAlbadaTVD);CHKERRQ(ierr);
596 ierr = PetscFunctionListAdd(&limiters,"koren" ,Limit_Koren);CHKERRQ(ierr);
597 ierr = PetscFunctionListAdd(&limiters,"korensym" ,Limit_KorenSym);CHKERRQ(ierr);
598 ierr = PetscFunctionListAdd(&limiters,"koren3" ,Limit_Koren3);CHKERRQ(ierr);
599 ierr = PetscFunctionListAdd(&limiters,"cada-torrilhon2" ,Limit_CadaTorrilhon2);CHKERRQ(ierr);
600 ierr = PetscFunctionListAdd(&limiters,"cada-torrilhon3-r0p1",Limit_CadaTorrilhon3R0p1);CHKERRQ(ierr);
601 ierr = PetscFunctionListAdd(&limiters,"cada-torrilhon3-r1" ,Limit_CadaTorrilhon3R1);CHKERRQ(ierr);
602 ierr = PetscFunctionListAdd(&limiters,"cada-torrilhon3-r10" ,Limit_CadaTorrilhon3R10);CHKERRQ(ierr);
603 ierr = PetscFunctionListAdd(&limiters,"cada-torrilhon3-r100",Limit_CadaTorrilhon3R100);CHKERRQ(ierr);
604
605 /* Register physical models to be available on the command line */
606 ierr = PetscFunctionListAdd(&physics,"advect" ,PhysicsCreate_Advect);CHKERRQ(ierr);
607
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 ierr = PetscOptionsBegin(comm,NULL,"Finite Volume solver options","");CHKERRQ(ierr);
614 ierr = PetscOptionsReal("-xmin","X min","",ctx.xmin,&ctx.xmin,NULL);CHKERRQ(ierr);
615 ierr = PetscOptionsReal("-xmax","X max","",ctx.xmax,&ctx.xmax,NULL);CHKERRQ(ierr);
616 ierr = PetscOptionsFList("-limit","Name of flux limiter to use","",limiters,lname,lname,sizeof(lname),NULL);CHKERRQ(ierr);
617 ierr = PetscOptionsInt("-draw","Draw solution vector, bitwise OR of (1=initial,2=final,4=final error)","",draw,&draw,NULL);CHKERRQ(ierr);
618 ierr = PetscOptionsString("-view_final","Write final solution in ASCII MATLAB format to given file name","",final_fname,final_fname,sizeof(final_fname),&view_final);CHKERRQ(ierr);
619 ierr = PetscOptionsInt("-initial","Initial condition (depends on the physics)","",ctx.initial,&ctx.initial,NULL);CHKERRQ(ierr);
620 ierr = PetscOptionsBool("-simulation","Compare errors with reference solution","",ctx.simulation,&ctx.simulation,NULL);CHKERRQ(ierr);
621 ierr = PetscOptionsReal("-cfl","CFL number to time step at","",ctx.cfl,&ctx.cfl,NULL);CHKERRQ(ierr);
622 ierr = PetscOptionsEnum("-bc_type","Boundary condition","",FVBCTypes,(PetscEnum)ctx.bctype,(PetscEnum*)&ctx.bctype,NULL);CHKERRQ(ierr);
623 ierr = PetscOptionsBool("-recursive_split","Split the domain recursively","",ctx.recursive,&ctx.recursive,NULL);CHKERRQ(ierr);
624 ierr = PetscOptionsEnd();CHKERRQ(ierr);
625
626 /* Choose the limiter from the list of registered limiters */
627 ierr = PetscFunctionListFind(limiters,lname,&ctx.limit);CHKERRQ(ierr);
628 if (!ctx.limit) SETERRQ1(PETSC_COMM_SELF,1,"Limiter '%s' not found",lname);
629
630 /* Choose the physics from the list of registered models */
631 {
632 PetscErrorCode (*r)(FVCtx*);
633 ierr = PetscFunctionListFind(physics,physname,&r);CHKERRQ(ierr);
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 ierr = (*r)(&ctx);CHKERRQ(ierr);
637 }
638
639 /* Create a DMDA to manage the parallel grid */
640 ierr = DMDACreate1d(comm,DM_BOUNDARY_PERIODIC,50,ctx.physics.dof,2,NULL,&da);CHKERRQ(ierr);
641 ierr = DMSetFromOptions(da);CHKERRQ(ierr);
642 ierr = DMSetUp(da);CHKERRQ(ierr);
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 ierr = DMDASetFieldName(da,i,ctx.physics.fieldname[i]);CHKERRQ(ierr);
647 }
648 ierr = DMDAGetInfo(da,0, &Mx,0,0, 0,0,0, &dof,0,0,0,0,0);CHKERRQ(ierr);
649 ierr = DMDAGetCorners(da,&xs,0,0,&xm,0,0);CHKERRQ(ierr);
650
651 /* Set coordinates of cell centers */
652 ierr = DMDASetUniformCoordinates(da,ctx.xmin+0.5*(ctx.xmax-ctx.xmin)/Mx,ctx.xmax+0.5*(ctx.xmax-ctx.xmin)/Mx,0,0,0,0);CHKERRQ(ierr);
653
654 /* Allocate work space for the Finite Volume solver (so it doesn't have to be reallocated on each function evaluation) */
655 ierr = PetscMalloc4(dof*dof,&ctx.R,dof*dof,&ctx.Rinv,2*dof,&ctx.cjmpLR,1*dof,&ctx.cslope);CHKERRQ(ierr);
656 ierr = PetscMalloc3(2*dof,&ctx.uLR,dof,&ctx.flux,dof,&ctx.speeds);CHKERRQ(ierr);
657
658 /* Create a vector to store the solution and to save the initial state */
659 ierr = DMCreateGlobalVector(da,&X);CHKERRQ(ierr);
660 ierr = VecDuplicate(X,&X0);CHKERRQ(ierr);
661 ierr = VecDuplicate(X,&R);CHKERRQ(ierr);
662
663 /*-------------------------------- create index for slow parts and fast parts ----------------------------------------*/
664 ierr = PetscMalloc1(xm*dof,&index_slow);CHKERRQ(ierr);
665 ierr = PetscMalloc1(xm*dof,&index_fast);CHKERRQ(ierr);
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 ierr = ISCreateGeneral(PETSC_COMM_WORLD,islow,index_slow,PETSC_COPY_VALUES,&ctx.iss);CHKERRQ(ierr);
673 ierr = ISCreateGeneral(PETSC_COMM_WORLD,ifast,index_fast,PETSC_COPY_VALUES,&ctx.isf);CHKERRQ(ierr);
674
675 /* Create a time-stepping object */
676 ierr = TSCreate(comm,&ts);CHKERRQ(ierr);
677 ierr = TSSetDM(ts,da);CHKERRQ(ierr);
678 ierr = TSSetRHSFunction(ts,R,FVRHSFunction,&ctx);CHKERRQ(ierr);
679 ierr = TSRHSSplitSetIS(ts,"slow",ctx.iss);CHKERRQ(ierr);
680 ierr = TSRHSSplitSetIS(ts,"fast",ctx.isf);CHKERRQ(ierr);
681 ierr = TSRHSSplitSetRHSFunction(ts,"slow",NULL,FVRHSFunctionslow,&ctx);CHKERRQ(ierr);
682 ierr = TSRHSSplitSetRHSFunction(ts,"fast",NULL,FVRHSFunctionfast,&ctx);CHKERRQ(ierr);
683
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 ierr = ISCreateGeneral(PETSC_COMM_WORLD,islow,index_slow,PETSC_COPY_VALUES,&ctx.iss2);CHKERRQ(ierr);
696 ierr = ISCreateGeneral(PETSC_COMM_WORLD,ifast,index_fast,PETSC_COPY_VALUES,&ctx.isf2);CHKERRQ(ierr);
697
698 ierr = TSRHSSplitGetSubTS(ts,"fast",&subts);CHKERRQ(ierr);
699 ierr = TSRHSSplitSetIS(subts,"slow",ctx.iss2);CHKERRQ(ierr);
700 ierr = TSRHSSplitSetIS(subts,"fast",ctx.isf2);CHKERRQ(ierr);
701 ierr = TSRHSSplitSetRHSFunction(subts,"slow",NULL,FVRHSFunctionslow2,&ctx);CHKERRQ(ierr);
702 ierr = TSRHSSplitSetRHSFunction(subts,"fast",NULL,FVRHSFunctionfast2,&ctx);CHKERRQ(ierr);
703 }
704
705 /*ierr = TSSetType(ts,TSSSP);CHKERRQ(ierr);*/
706 ierr = TSSetType(ts,TSMPRK);CHKERRQ(ierr);
707 ierr = TSSetMaxTime(ts,10);CHKERRQ(ierr);
708 ierr = TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);CHKERRQ(ierr);
709
710 /* Compute initial conditions and starting time step */
711 ierr = FVSample(&ctx,da,0,X0);CHKERRQ(ierr);
712 ierr = FVRHSFunction(ts,0,X0,X,(void*)&ctx);CHKERRQ(ierr); /* Initial function evaluation, only used to determine max speed */
713 ierr = VecCopy(X0,X);CHKERRQ(ierr); /* The function value was not used so we set X=X0 again */
714 ierr = TSSetTimeStep(ts,ctx.cfl/ctx.cfl_idt);CHKERRQ(ierr);
715 ierr = TSSetFromOptions(ts);CHKERRQ(ierr); /* Take runtime options */
716 ierr = SolutionStatsView(da,X,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
717 {
718 PetscInt steps;
719 PetscScalar mass_initial,mass_final,mass_difference;
720
721 ierr = TSSolve(ts,X);CHKERRQ(ierr);
722 ierr = TSGetSolveTime(ts,&ptime);CHKERRQ(ierr);
723 ierr = TSGetStepNumber(ts,&steps);CHKERRQ(ierr);
724 ierr = PetscPrintf(comm,"Final time %g, steps %D\n",(double)ptime,steps);CHKERRQ(ierr);
725 /* calculate the total mass at initial time and final time */
726 mass_initial = 0.0;
727 mass_final = 0.0;
728 ierr = VecSum(X0,&mass_initial);CHKERRQ(ierr);
729 ierr = VecSum(X,&mass_final);CHKERRQ(ierr);
730 mass_difference = (ctx.xmax-ctx.xmin)/(PetscScalar)Mx*(mass_final - mass_initial);
731 ierr = PetscPrintf(comm,"Mass difference %g\n",(double)mass_difference);CHKERRQ(ierr);
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 ierr = PetscOptionsGetString(NULL,NULL,"-f",filename,sizeof(filename),&flg);CHKERRQ(ierr);
739 if (!flg) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_USER,"Must indicate binary file with the -f option");
740 ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_READ,&fd);CHKERRQ(ierr);
741 ierr = VecDuplicate(X0,&XR);CHKERRQ(ierr);
742 ierr = VecLoad(XR,fd);CHKERRQ(ierr);
743 ierr = PetscViewerDestroy(&fd);CHKERRQ(ierr);
744 ierr = VecAYPX(XR,-1,X);CHKERRQ(ierr);
745 ierr = VecNorm(XR,NORM_1,&nrm1);CHKERRQ(ierr);
746 ierr = VecNorm(XR,NORM_INFINITY,&nrmsup);CHKERRQ(ierr);
747 nrm1 /= Mx;
748 ierr = PetscPrintf(comm,"Error ||x-x_e||_1 %g ||x-x_e||_sup %g\n",(double)nrm1,(double)nrmsup);CHKERRQ(ierr);
749 ierr = VecDestroy(&XR);CHKERRQ(ierr);
750 }
751 }
752
753 ierr = SolutionStatsView(da,X,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
754 if (draw & 0x1) {ierr = VecView(X0,PETSC_VIEWER_DRAW_WORLD);CHKERRQ(ierr);}
755 if (draw & 0x2) {ierr = VecView(X,PETSC_VIEWER_DRAW_WORLD);CHKERRQ(ierr);}
756 if (draw & 0x4) {
757 Vec Y;
758 ierr = VecDuplicate(X,&Y);CHKERRQ(ierr);
759 ierr = FVSample(&ctx,da,ptime,Y);CHKERRQ(ierr);
760 ierr = VecAYPX(Y,-1,X);CHKERRQ(ierr);
761 ierr = VecView(Y,PETSC_VIEWER_DRAW_WORLD);CHKERRQ(ierr);
762 ierr = VecDestroy(&Y);CHKERRQ(ierr);
763 }
764
765 if (view_final) {
766 PetscViewer viewer;
767 ierr = PetscViewerASCIIOpen(PETSC_COMM_WORLD,final_fname,&viewer);CHKERRQ(ierr);
768 ierr = PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_MATLAB);CHKERRQ(ierr);
769 ierr = VecView(X,viewer);CHKERRQ(ierr);
770 ierr = PetscViewerPopFormat(viewer);CHKERRQ(ierr);
771 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);
772 }
773
774 /* Clean up */
775 ierr = (*ctx.physics.destroy)(ctx.physics.user);CHKERRQ(ierr);
776 for (i=0; i<ctx.physics.dof; i++) {ierr = PetscFree(ctx.physics.fieldname[i]);CHKERRQ(ierr);}
777 ierr = PetscFree4(ctx.R,ctx.Rinv,ctx.cjmpLR,ctx.cslope);CHKERRQ(ierr);
778 ierr = PetscFree3(ctx.uLR,ctx.flux,ctx.speeds);CHKERRQ(ierr);
779 ierr = ISDestroy(&ctx.iss);CHKERRQ(ierr);
780 ierr = ISDestroy(&ctx.isf);CHKERRQ(ierr);
781 ierr = ISDestroy(&ctx.iss2);CHKERRQ(ierr);
782 ierr = ISDestroy(&ctx.isf2);CHKERRQ(ierr);
783 ierr = VecDestroy(&X);CHKERRQ(ierr);
784 ierr = VecDestroy(&X0);CHKERRQ(ierr);
785 ierr = VecDestroy(&R);CHKERRQ(ierr);
786 ierr = DMDestroy(&da);CHKERRQ(ierr);
787 ierr = TSDestroy(&ts);CHKERRQ(ierr);
788 ierr = PetscFree(index_slow);CHKERRQ(ierr);
789 ierr = PetscFree(index_fast);CHKERRQ(ierr);
790 ierr = PetscFunctionListDestroy(&limiters);CHKERRQ(ierr);
791 ierr = PetscFunctionListDestroy(&physics);CHKERRQ(ierr);
792 ierr = PetscFinalize();
793 return ierr;
794 }
795
796 /*TEST
797 build:
798 requires: !complex
799 depends: finitevolume1d.c
800
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
803
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
808
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
812
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
817
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
821
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*/
826