1 /*
2 Note:
3 -hratio is the ratio between mesh size of corse grids and fine grids
4 -ts_rk_dtratio is the ratio between the large stepsize and the small stepsize
5 */
6
7 static const char help[] = "1D periodic Finite Volume solver in slope-limiter form with semidiscrete time stepping.\n"
8 " advection - Constant coefficient scalar advection\n"
9 " u_t + (a*u)_x = 0\n"
10 " for this toy problem, we choose different meshsizes for different sub-domains (slow-fast-slow), say\n"
11 " hxs = hratio*hxf \n"
12 " where hxs and hxf are the grid spacings for coarse and fine grids respectively.\n"
13 " exact - Exact Riemann solver which usually needs to perform a Newton iteration to connect\n"
14 " the states across shocks and rarefactions\n"
15 " simulation - use reference solution which is generated by smaller time step size to be true solution,\n"
16 " also the reference solution should be generated by user and stored in a binary file.\n"
17 " characteristic - Limit the characteristic variables, this is usually preferred (default)\n"
18 "Several initial conditions can be chosen with -initial N\n\n"
19 "The problem size should be set with -da_grid_x M\n\n";
20
21 #include <petscts.h>
22 #include <petscdm.h>
23 #include <petscdmda.h>
24 #include <petscdraw.h>
25 #include "finitevolume1d.h"
26
RangeMod(PetscReal a,PetscReal xmin,PetscReal xmax)27 PETSC_STATIC_INLINE PetscReal RangeMod(PetscReal a,PetscReal xmin,PetscReal xmax) { PetscReal range = xmax-xmin; return xmin +PetscFmodReal(range+PetscFmodReal(a,range),range); }
28
29 /* --------------------------------- Advection ----------------------------------- */
30 typedef struct {
31 PetscReal a; /* advective velocity */
32 } AdvectCtx;
33
PhysicsRiemann_Advect(void * vctx,PetscInt m,const PetscScalar * uL,const PetscScalar * uR,PetscScalar * flux,PetscReal * maxspeed)34 static PetscErrorCode PhysicsRiemann_Advect(void *vctx,PetscInt m,const PetscScalar *uL,const PetscScalar *uR,PetscScalar *flux,PetscReal *maxspeed)
35 {
36 AdvectCtx *ctx = (AdvectCtx*)vctx;
37 PetscReal speed;
38
39 PetscFunctionBeginUser;
40 speed = ctx->a;
41 flux[0] = PetscMax(0,speed)*uL[0] + PetscMin(0,speed)*uR[0];
42 *maxspeed = speed;
43 PetscFunctionReturn(0);
44 }
45
PhysicsCharacteristic_Advect(void * vctx,PetscInt m,const PetscScalar * u,PetscScalar * X,PetscScalar * Xi,PetscReal * speeds)46 static PetscErrorCode PhysicsCharacteristic_Advect(void *vctx,PetscInt m,const PetscScalar *u,PetscScalar *X,PetscScalar *Xi,PetscReal *speeds)
47 {
48 AdvectCtx *ctx = (AdvectCtx*)vctx;
49
50 PetscFunctionBeginUser;
51 X[0] = 1.;
52 Xi[0] = 1.;
53 speeds[0] = ctx->a;
54 PetscFunctionReturn(0);
55 }
56
PhysicsSample_Advect(void * vctx,PetscInt initial,FVBCType bctype,PetscReal xmin,PetscReal xmax,PetscReal t,PetscReal x,PetscReal * u)57 static PetscErrorCode PhysicsSample_Advect(void *vctx,PetscInt initial,FVBCType bctype,PetscReal xmin,PetscReal xmax,PetscReal t,PetscReal x,PetscReal *u)
58 {
59 AdvectCtx *ctx = (AdvectCtx*)vctx;
60 PetscReal a = ctx->a,x0;
61
62 PetscFunctionBeginUser;
63 switch (bctype) {
64 case FVBC_OUTFLOW: x0 = x-a*t; break;
65 case FVBC_PERIODIC: x0 = RangeMod(x-a*t,xmin,xmax); break;
66 default: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"unknown BCType");
67 }
68 switch (initial) {
69 case 0: u[0] = (x0 < 0) ? 1 : -1; break;
70 case 1: u[0] = (x0 < 0) ? -1 : 1; break;
71 case 2: u[0] = (0 < x0 && x0 < 1) ? 1 : 0; break;
72 case 3: u[0] = PetscSinReal(2*PETSC_PI*x0); break;
73 case 4: u[0] = PetscAbs(x0); break;
74 case 5: u[0] = (x0 < 0 || x0 > 0.5) ? 0 : PetscSqr(PetscSinReal(2*PETSC_PI*x0)); break;
75 case 6: u[0] = (x0 < 0) ? 0 : ((x0 < 1) ? x0 : ((x0 < 2) ? 2-x0 : 0)); break;
76 case 7: u[0] = PetscPowReal(PetscSinReal(PETSC_PI*x0),10.0);break;
77 default: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"unknown initial condition");
78 }
79 PetscFunctionReturn(0);
80 }
81
PhysicsCreate_Advect(FVCtx * ctx)82 static PetscErrorCode PhysicsCreate_Advect(FVCtx *ctx)
83 {
84 PetscErrorCode ierr;
85 AdvectCtx *user;
86
87 PetscFunctionBeginUser;
88 ierr = PetscNew(&user);CHKERRQ(ierr);
89 ctx->physics2.sample2 = PhysicsSample_Advect;
90 ctx->physics2.riemann2 = PhysicsRiemann_Advect;
91 ctx->physics2.characteristic2 = PhysicsCharacteristic_Advect;
92 ctx->physics2.destroy = PhysicsDestroy_SimpleFree;
93 ctx->physics2.user = user;
94 ctx->physics2.dof = 1;
95 ierr = PetscStrallocpy("u",&ctx->physics2.fieldname[0]);CHKERRQ(ierr);
96 user->a = 1;
97 ierr = PetscOptionsBegin(ctx->comm,ctx->prefix,"Options for advection","");CHKERRQ(ierr);
98 {
99 ierr = PetscOptionsReal("-physics_advect_a","Speed","",user->a,&user->a,NULL);CHKERRQ(ierr);
100 }
101 ierr = PetscOptionsEnd();CHKERRQ(ierr);
102 PetscFunctionReturn(0);
103 }
104
FVSample_2WaySplit(FVCtx * ctx,DM da,PetscReal time,Vec U)105 PetscErrorCode FVSample_2WaySplit(FVCtx *ctx,DM da,PetscReal time,Vec U)
106 {
107 PetscErrorCode ierr;
108 PetscScalar *u,*uj,xj,xi;
109 PetscInt i,j,k,dof,xs,xm,Mx;
110 const PetscInt N = 200;
111 PetscReal hs,hf;
112
113 PetscFunctionBeginUser;
114 if (!ctx->physics2.sample2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Physics has not provided a sampling function");
115 ierr = DMDAGetInfo(da,0, &Mx,0,0, 0,0,0, &dof,0,0,0,0,0);CHKERRQ(ierr);
116 ierr = DMDAGetCorners(da,&xs,0,0,&xm,0,0);CHKERRQ(ierr);
117 ierr = DMDAVecGetArray(da,U,&u);CHKERRQ(ierr);
118 ierr = PetscMalloc1(dof,&uj);CHKERRQ(ierr);
119 hs = (ctx->xmax-ctx->xmin)*3.0/8.0/ctx->sf;
120 hf = (ctx->xmax-ctx->xmin)/4.0/(ctx->fs-ctx->sf);
121 for (i=xs; i<xs+xm; i++) {
122 if (i < ctx->sf) {
123 xi = ctx->xmin+0.5*hs+i*hs;
124 /* Integrate over cell i using trapezoid rule with N points. */
125 for (k=0; k<dof; k++) u[i*dof+k] = 0;
126 for (j=0; j<N+1; j++) {
127 xj = xi+hs*(j-N/2)/(PetscReal)N;
128 ierr = (*ctx->physics2.sample2)(ctx->physics2.user,ctx->initial,ctx->bctype,ctx->xmin,ctx->xmax,time,xj,uj);CHKERRQ(ierr);
129 for (k=0; k<dof; k++) u[i*dof+k] += ((j==0 || j==N) ? 0.5 : 1.0)*uj[k]/N;
130 }
131 } else if (i < ctx->fs) {
132 xi = ctx->xmin+ctx->sf*hs+0.5*hf+(i-ctx->sf)*hf;
133 /* Integrate over cell i using trapezoid rule with N points. */
134 for (k=0; k<dof; k++) u[i*dof+k] = 0;
135 for (j=0; j<N+1; j++) {
136 xj = xi+hf*(j-N/2)/(PetscReal)N;
137 ierr = (*ctx->physics2.sample2)(ctx->physics2.user,ctx->initial,ctx->bctype,ctx->xmin,ctx->xmax,time,xj,uj);CHKERRQ(ierr);
138 for (k=0; k<dof; k++) u[i*dof+k] += ((j==0 || j==N) ? 0.5 : 1.0)*uj[k]/N;
139 }
140 } else {
141 xi = ctx->xmin+ctx->sf*hs+(ctx->fs-ctx->sf)*hf+0.5*hs+(i-ctx->fs)*hs;
142 /* Integrate over cell i using trapezoid rule with N points. */
143 for (k=0; k<dof; k++) u[i*dof+k] = 0;
144 for (j=0; j<N+1; j++) {
145 xj = xi+hs*(j-N/2)/(PetscReal)N;
146 ierr = (*ctx->physics2.sample2)(ctx->physics2.user,ctx->initial,ctx->bctype,ctx->xmin,ctx->xmax,time,xj,uj);CHKERRQ(ierr);
147 for (k=0; k<dof; k++) u[i*dof+k] += ((j==0 || j==N) ? 0.5 : 1.0)*uj[k]/N;
148 }
149 }
150 }
151 ierr = DMDAVecRestoreArray(da,U,&u);CHKERRQ(ierr);
152 ierr = PetscFree(uj);CHKERRQ(ierr);
153 PetscFunctionReturn(0);
154 }
155
SolutionErrorNorms_2WaySplit(FVCtx * ctx,DM da,PetscReal t,Vec X,PetscReal * nrm1)156 static PetscErrorCode SolutionErrorNorms_2WaySplit(FVCtx *ctx,DM da,PetscReal t,Vec X,PetscReal *nrm1)
157 {
158 PetscErrorCode ierr;
159 Vec Y;
160 PetscInt i,Mx;
161 const PetscScalar *ptr_X,*ptr_Y;
162 PetscReal hs,hf;
163
164 PetscFunctionBeginUser;
165 ierr = VecGetSize(X,&Mx);CHKERRQ(ierr);
166 ierr = VecDuplicate(X,&Y);CHKERRQ(ierr);
167 ierr = FVSample_2WaySplit(ctx,da,t,Y);CHKERRQ(ierr);
168 hs = (ctx->xmax-ctx->xmin)*3.0/8.0/ctx->sf;
169 hf = (ctx->xmax-ctx->xmin)/4.0/(ctx->fs-ctx->sf);
170 ierr = VecGetArrayRead(X,&ptr_X);CHKERRQ(ierr);
171 ierr = VecGetArrayRead(Y,&ptr_Y);CHKERRQ(ierr);
172 for (i=0; i<Mx; i++) {
173 if (i < ctx->sf || i > ctx->fs -1) *nrm1 += hs*PetscAbs(ptr_X[i]-ptr_Y[i]);
174 else *nrm1 += hf*PetscAbs(ptr_X[i]-ptr_Y[i]);
175 }
176 ierr = VecRestoreArrayRead(X,&ptr_X);CHKERRQ(ierr);
177 ierr = VecRestoreArrayRead(Y,&ptr_Y);CHKERRQ(ierr);
178 ierr = VecDestroy(&Y);CHKERRQ(ierr);
179 PetscFunctionReturn(0);
180 }
181
FVRHSFunction_2WaySplit(TS ts,PetscReal time,Vec X,Vec F,void * vctx)182 PetscErrorCode FVRHSFunction_2WaySplit(TS ts,PetscReal time,Vec X,Vec F,void *vctx)
183 {
184 FVCtx *ctx = (FVCtx*)vctx;
185 PetscErrorCode ierr;
186 PetscInt i,j,k,Mx,dof,xs,xm,sf = ctx->sf,fs = ctx->fs;
187 PetscReal hxf,hxs,cfl_idt = 0;
188 PetscScalar *x,*f,*slope;
189 Vec Xloc;
190 DM da;
191
192 PetscFunctionBeginUser;
193 ierr = TSGetDM(ts,&da);CHKERRQ(ierr);
194 ierr = DMGetLocalVector(da,&Xloc);CHKERRQ(ierr); /* Xloc contains ghost points */
195 ierr = DMDAGetInfo(da,0, &Mx,0,0, 0,0,0, &dof,0,0,0,0,0);CHKERRQ(ierr); /* Mx is the number of center points */
196 hxs = (ctx->xmax-ctx->xmin)*3.0/8.0/ctx->sf;
197 hxf = (ctx->xmax-ctx->xmin)/4.0/(ctx->fs-ctx->sf);
198 ierr = DMGlobalToLocalBegin(da,X,INSERT_VALUES,Xloc);CHKERRQ(ierr); /* X is solution vector which does not contain ghost points */
199 ierr = DMGlobalToLocalEnd(da,X,INSERT_VALUES,Xloc);CHKERRQ(ierr);
200
201 ierr = VecZeroEntries(F);CHKERRQ(ierr); /* F is the right hand side function corresponds to center points */
202
203 ierr = DMDAVecGetArray(da,Xloc,&x);CHKERRQ(ierr);
204 ierr = DMDAVecGetArray(da,F,&f);CHKERRQ(ierr);
205 ierr = DMDAGetArray(da,PETSC_TRUE,&slope);CHKERRQ(ierr); /* contains ghost points */
206
207 ierr = DMDAGetCorners(da,&xs,0,0,&xm,0,0);CHKERRQ(ierr);
208
209 if (ctx->bctype == FVBC_OUTFLOW) {
210 for (i=xs-2; i<0; i++) {
211 for (j=0; j<dof; j++) x[i*dof+j] = x[j];
212 }
213 for (i=Mx; i<xs+xm+2; i++) {
214 for (j=0; j<dof; j++) x[i*dof+j] = x[(xs+xm-1)*dof+j];
215 }
216 }
217 for (i=xs-1; i<xs+xm+1; i++) {
218 struct _LimitInfo info;
219 PetscScalar *cjmpL,*cjmpR;
220 /* Determine the right eigenvectors R, where A = R \Lambda R^{-1} */
221 ierr = (*ctx->physics2.characteristic2)(ctx->physics2.user,dof,&x[i*dof],ctx->R,ctx->Rinv,ctx->speeds);CHKERRQ(ierr);
222 /* Evaluate jumps across interfaces (i-1, i) and (i, i+1), put in characteristic basis */
223 ierr = PetscArrayzero(ctx->cjmpLR,2*dof);CHKERRQ(ierr);
224 cjmpL = &ctx->cjmpLR[0];
225 cjmpR = &ctx->cjmpLR[dof];
226 for (j=0; j<dof; j++) {
227 PetscScalar jmpL,jmpR;
228 jmpL = x[(i+0)*dof+j]-x[(i-1)*dof+j];
229 jmpR = x[(i+1)*dof+j]-x[(i+0)*dof+j];
230 for (k=0; k<dof; k++) {
231 cjmpL[k] += ctx->Rinv[k+j*dof]*jmpL;
232 cjmpR[k] += ctx->Rinv[k+j*dof]*jmpR;
233 }
234 }
235 /* Apply limiter to the left and right characteristic jumps */
236 info.m = dof;
237 info.hxs = hxs;
238 info.hxf = hxf;
239 (*ctx->limit2)(&info,cjmpL,cjmpR,ctx->sf,ctx->fs,i,ctx->cslope);
240 for (j=0; j<dof; j++) {
241 PetscScalar tmp = 0;
242 for (k=0; k<dof; k++) tmp += ctx->R[j+k*dof]*ctx->cslope[k];
243 slope[i*dof+j] = tmp;
244 }
245 }
246
247 for (i=xs; i<xs+xm+1; i++) {
248 PetscReal maxspeed;
249 PetscScalar *uL,*uR;
250 uL = &ctx->uLR[0];
251 uR = &ctx->uLR[dof];
252 if (i < sf) { /* slow region */
253 for (j=0; j<dof; j++) {
254 uL[j] = x[(i-1)*dof+j]+slope[(i-1)*dof+j]*hxs/2;
255 uR[j] = x[(i-0)*dof+j]-slope[(i-0)*dof+j]*hxs/2;
256 }
257 ierr = (*ctx->physics2.riemann2)(ctx->physics2.user,dof,uL,uR,ctx->flux,&maxspeed);CHKERRQ(ierr);
258 if (i > xs) {
259 for (j=0; j<dof; j++) f[(i-1)*dof+j] -= ctx->flux[j]/hxs;
260 }
261 if (i < xs+xm) {
262 for (j=0; j<dof; j++) f[i*dof+j] += ctx->flux[j]/hxs;
263 }
264 } else if (i == sf) { /* interface between the slow region and the fast region */
265 for (j=0; j<dof; j++) {
266 uL[j] = x[(i-1)*dof+j]+slope[(i-1)*dof+j]*hxs/2;
267 uR[j] = x[(i-0)*dof+j]-slope[(i-0)*dof+j]*hxf/2;
268 }
269 ierr = (*ctx->physics2.riemann2)(ctx->physics2.user,dof,uL,uR,ctx->flux,&maxspeed);CHKERRQ(ierr);
270 if (i > xs) {
271 for (j=0; j<dof; j++) f[(i-1)*dof+j] -= ctx->flux[j]/hxs;
272 }
273 if (i < xs+xm) {
274 for (j=0; j<dof; j++) f[i*dof+j] += ctx->flux[j]/hxf;
275 }
276 } else if (i > sf && i < fs) { /* fast region */
277 for (j=0; j<dof; j++) {
278 uL[j] = x[(i-1)*dof+j]+slope[(i-1)*dof+j]*hxf/2;
279 uR[j] = x[(i-0)*dof+j]-slope[(i-0)*dof+j]*hxf/2;
280 }
281 ierr = (*ctx->physics2.riemann2)(ctx->physics2.user,dof,uL,uR,ctx->flux,&maxspeed);CHKERRQ(ierr);
282 if (i > xs) {
283 for (j=0; j<dof; j++) f[(i-1)*dof+j] -= ctx->flux[j]/hxf;
284 }
285 if (i < xs+xm) {
286 for (j=0; j<dof; j++) f[i*dof+j] += ctx->flux[j]/hxf;
287 }
288 } else if (i == fs) { /* interface between the fast region and the slow region */
289 for (j=0; j<dof; j++) {
290 uL[j] = x[(i-1)*dof+j]+slope[(i-1)*dof+j]*hxf/2;
291 uR[j] = x[(i-0)*dof+j]-slope[(i-0)*dof+j]*hxs/2;
292 }
293 ierr = (*ctx->physics2.riemann2)(ctx->physics2.user,dof,uL,uR,ctx->flux,&maxspeed);CHKERRQ(ierr);
294 if (i > xs) {
295 for (j=0; j<dof; j++) f[(i-1)*dof+j] -= ctx->flux[j]/hxf;
296 }
297 if (i < xs+xm) {
298 for (j=0; j<dof; j++) f[i*dof+j] += ctx->flux[j]/hxs;
299 }
300 } else { /* slow region */
301 for (j=0; j<dof; j++) {
302 uL[j] = x[(i-1)*dof+j]+slope[(i-1)*dof+j]*hxs/2;
303 uR[j] = x[(i-0)*dof+j]-slope[(i-0)*dof+j]*hxs/2;
304 }
305 ierr = (*ctx->physics2.riemann2)(ctx->physics2.user,dof,uL,uR,ctx->flux,&maxspeed);CHKERRQ(ierr);
306 cfl_idt = PetscMax(cfl_idt,PetscAbsScalar(maxspeed/hxs)); /* Max allowable value of 1/Delta t */
307 if (i > xs) {
308 for (j=0; j<dof; j++) f[(i-1)*dof+j] -= ctx->flux[j]/hxs;
309 }
310 if (i < xs+xm) {
311 for (j=0; j<dof; j++) f[i*dof+j] += ctx->flux[j]/hxs;
312 }
313 }
314 }
315 ierr = DMDAVecRestoreArray(da,Xloc,&x);CHKERRQ(ierr);
316 ierr = DMDAVecRestoreArray(da,F,&f);CHKERRQ(ierr);
317 ierr = DMDARestoreArray(da,PETSC_TRUE,&slope);CHKERRQ(ierr);
318 ierr = DMRestoreLocalVector(da,&Xloc);CHKERRQ(ierr);
319 ierr = MPI_Allreduce(&cfl_idt,&ctx->cfl_idt,1,MPIU_SCALAR,MPIU_MAX,PetscObjectComm((PetscObject)da));CHKERRQ(ierr);
320 if (0) {
321 /* We need to a way to inform the TS of a CFL constraint, this is a debugging fragment */
322 PetscReal dt,tnow;
323 ierr = TSGetTimeStep(ts,&dt);CHKERRQ(ierr);
324 ierr = TSGetTime(ts,&tnow);CHKERRQ(ierr);
325 if (dt > 0.5/ctx->cfl_idt) {
326 if (1) {
327 ierr = PetscPrintf(ctx->comm,"Stability constraint exceeded at t=%g, dt %g > %g\n",(double)tnow,(double)dt,(double)(0.5/ctx->cfl_idt));CHKERRQ(ierr);
328 } else SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Stability constraint exceeded, %g > %g",(double)dt,(double)(ctx->cfl/ctx->cfl_idt));
329 }
330 }
331 PetscFunctionReturn(0);
332 }
333
334 /* --------------------------------- Finite Volume Solver for slow components ----------------------------------- */
FVRHSFunctionslow_2WaySplit(TS ts,PetscReal time,Vec X,Vec F,void * vctx)335 PetscErrorCode FVRHSFunctionslow_2WaySplit(TS ts,PetscReal time,Vec X,Vec F,void *vctx)
336 {
337 FVCtx *ctx = (FVCtx*)vctx;
338 PetscErrorCode ierr;
339 PetscInt i,j,k,Mx,dof,xs,xm,islow = 0,sf = ctx->sf,fs = ctx->fs,lsbwidth = ctx->lsbwidth,rsbwidth = ctx->rsbwidth;
340 PetscReal hxs,hxf,cfl_idt = 0;
341 PetscScalar *x,*f,*slope;
342 Vec Xloc;
343 DM da;
344
345 PetscFunctionBeginUser;
346 ierr = TSGetDM(ts,&da);CHKERRQ(ierr);
347 ierr = DMGetLocalVector(da,&Xloc);CHKERRQ(ierr);
348 ierr = DMDAGetInfo(da,0, &Mx,0,0, 0,0,0, &dof,0,0,0,0,0);CHKERRQ(ierr);
349 hxs = (ctx->xmax-ctx->xmin)*3.0/8.0/ctx->sf;
350 hxf = (ctx->xmax-ctx->xmin)/4.0/(ctx->fs-ctx->sf);
351 ierr = DMGlobalToLocalBegin(da,X,INSERT_VALUES,Xloc);CHKERRQ(ierr);
352 ierr = DMGlobalToLocalEnd(da,X,INSERT_VALUES,Xloc);CHKERRQ(ierr);
353 ierr = VecZeroEntries(F);CHKERRQ(ierr);
354 ierr = DMDAVecGetArray(da,Xloc,&x);CHKERRQ(ierr);
355 ierr = VecGetArray(F,&f);CHKERRQ(ierr);
356 ierr = DMDAGetArray(da,PETSC_TRUE,&slope);CHKERRQ(ierr);
357 ierr = DMDAGetCorners(da,&xs,0,0,&xm,0,0);CHKERRQ(ierr);
358
359 if (ctx->bctype == FVBC_OUTFLOW) {
360 for (i=xs-2; i<0; i++) {
361 for (j=0; j<dof; j++) x[i*dof+j] = x[j];
362 }
363 for (i=Mx; i<xs+xm+2; i++) {
364 for (j=0; j<dof; j++) x[i*dof+j] = x[(xs+xm-1)*dof+j];
365 }
366 }
367 for (i=xs-1; i<xs+xm+1; i++) {
368 struct _LimitInfo info;
369 PetscScalar *cjmpL,*cjmpR;
370 if (i < sf-lsbwidth+1 || i > fs+rsbwidth-2) { /* slow components and the first and last fast components */
371 /* Determine the right eigenvectors R, where A = R \Lambda R^{-1} */
372 ierr = (*ctx->physics2.characteristic2)(ctx->physics2.user,dof,&x[i*dof],ctx->R,ctx->Rinv,ctx->speeds);CHKERRQ(ierr);
373 /* Evaluate jumps across interfaces (i-1, i) and (i, i+1), put in characteristic basis */
374 ierr = PetscArrayzero(ctx->cjmpLR,2*dof);CHKERRQ(ierr);
375 cjmpL = &ctx->cjmpLR[0];
376 cjmpR = &ctx->cjmpLR[dof];
377 for (j=0; j<dof; j++) {
378 PetscScalar jmpL,jmpR;
379 jmpL = x[(i+0)*dof+j]-x[(i-1)*dof+j];
380 jmpR = x[(i+1)*dof+j]-x[(i+0)*dof+j];
381 for (k=0; k<dof; k++) {
382 cjmpL[k] += ctx->Rinv[k+j*dof]*jmpL;
383 cjmpR[k] += ctx->Rinv[k+j*dof]*jmpR;
384 }
385 }
386 /* Apply limiter to the left and right characteristic jumps */
387 info.m = dof;
388 info.hxs = hxs;
389 info.hxf = hxf;
390 (*ctx->limit2)(&info,cjmpL,cjmpR,ctx->sf,ctx->fs,i,ctx->cslope);
391 for (j=0; j<dof; j++) {
392 PetscScalar tmp = 0;
393 for (k=0; k<dof; k++) tmp += ctx->R[j+k*dof]*ctx->cslope[k];
394 slope[i*dof+j] = tmp;
395 }
396 }
397 }
398
399 for (i=xs; i<xs+xm+1; i++) {
400 PetscReal maxspeed;
401 PetscScalar *uL,*uR;
402 uL = &ctx->uLR[0];
403 uR = &ctx->uLR[dof];
404 if (i < sf-lsbwidth) { /* slow region */
405 for (j=0; j<dof; j++) {
406 uL[j] = x[(i-1)*dof+j]+slope[(i-1)*dof+j]*hxs/2;
407 uR[j] = x[(i-0)*dof+j]-slope[(i-0)*dof+j]*hxs/2;
408 }
409 ierr = (*ctx->physics2.riemann2)(ctx->physics2.user,dof,uL,uR,ctx->flux,&maxspeed);CHKERRQ(ierr);
410 cfl_idt = PetscMax(cfl_idt,PetscAbsScalar(maxspeed/hxs)); /* Max allowable value of 1/Delta t */
411 if (i > xs) {
412 for (j=0; j<dof; j++) f[(islow-1)*dof+j] -= ctx->flux[j]/hxs;
413 }
414 if (i < xs+xm) {
415 for (j=0; j<dof; j++) f[islow*dof+j] += ctx->flux[j]/hxs;
416 islow++;
417 }
418 }
419 if (i == sf-lsbwidth) { /* interface between the slow region and the fast region */
420 for (j=0; j<dof; j++) {
421 uL[j] = x[(i-1)*dof+j]+slope[(i-1)*dof+j]*hxs/2;
422 uR[j] = x[(i-0)*dof+j]-slope[(i-0)*dof+j]*hxs/2;
423 }
424 ierr = (*ctx->physics2.riemann2)(ctx->physics2.user,dof,uL,uR,ctx->flux,&maxspeed);CHKERRQ(ierr);
425 if (i > xs) {
426 for (j=0; j<dof; j++) f[(islow-1)*dof+j] -= ctx->flux[j]/hxs;
427 }
428 }
429 if (i == fs+rsbwidth) { /* slow region */
430 for (j=0; j<dof; j++) {
431 uL[j] = x[(i-1)*dof+j]+slope[(i-1)*dof+j]*hxs/2;
432 uR[j] = x[(i-0)*dof+j]-slope[(i-0)*dof+j]*hxs/2;
433 }
434 ierr = (*ctx->physics2.riemann2)(ctx->physics2.user,dof,uL,uR,ctx->flux,&maxspeed);CHKERRQ(ierr);
435 if (i < xs+xm) {
436 for (j=0; j<dof; j++) f[islow*dof+j] += ctx->flux[j]/hxs;
437 islow++;
438 }
439 }
440 if (i > fs+rsbwidth) { /* slow region */
441 for (j=0; j<dof; j++) {
442 uL[j] = x[(i-1)*dof+j]+slope[(i-1)*dof+j]*hxs/2;
443 uR[j] = x[(i-0)*dof+j]-slope[(i-0)*dof+j]*hxs/2;
444 }
445 ierr = (*ctx->physics2.riemann2)(ctx->physics2.user,dof,uL,uR,ctx->flux,&maxspeed);CHKERRQ(ierr);
446 if (i > xs) {
447 for (j=0; j<dof; j++) f[(islow-1)*dof+j] -= ctx->flux[j]/hxs;
448 }
449 if (i < xs+xm) {
450 for (j=0; j<dof; j++) f[islow*dof+j] += ctx->flux[j]/hxs;
451 islow++;
452 }
453 }
454 }
455 ierr = DMDAVecRestoreArray(da,Xloc,&x);CHKERRQ(ierr);
456 ierr = VecRestoreArray(F,&f);CHKERRQ(ierr);
457 ierr = DMDARestoreArray(da,PETSC_TRUE,&slope);CHKERRQ(ierr);
458 ierr = DMRestoreLocalVector(da,&Xloc);CHKERRQ(ierr);
459 ierr = MPI_Allreduce(&cfl_idt,&ctx->cfl_idt,1,MPIU_SCALAR,MPIU_MAX,PetscObjectComm((PetscObject)da));CHKERRQ(ierr);
460 PetscFunctionReturn(0);
461 }
462
FVRHSFunctionslowbuffer_2WaySplit(TS ts,PetscReal time,Vec X,Vec F,void * vctx)463 PetscErrorCode FVRHSFunctionslowbuffer_2WaySplit(TS ts,PetscReal time,Vec X,Vec F,void *vctx)
464 {
465 FVCtx *ctx = (FVCtx*)vctx;
466 PetscErrorCode ierr;
467 PetscInt i,j,k,Mx,dof,xs,xm,islow = 0,sf = ctx->sf,fs = ctx->fs,lsbwidth = ctx->lsbwidth,rsbwidth = ctx->rsbwidth;
468 PetscReal hxs,hxf;
469 PetscScalar *x,*f,*slope;
470 Vec Xloc;
471 DM da;
472
473 PetscFunctionBeginUser;
474 ierr = TSGetDM(ts,&da);CHKERRQ(ierr);
475 ierr = DMGetLocalVector(da,&Xloc);CHKERRQ(ierr);
476 ierr = DMDAGetInfo(da,0, &Mx,0,0, 0,0,0, &dof,0,0,0,0,0);CHKERRQ(ierr);
477 hxs = (ctx->xmax-ctx->xmin)*3.0/8.0/ctx->sf;
478 hxf = (ctx->xmax-ctx->xmin)/4.0/(ctx->fs-ctx->sf);
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
487 if (ctx->bctype == FVBC_OUTFLOW) {
488 for (i=xs-2; i<0; i++) {
489 for (j=0; j<dof; j++) x[i*dof+j] = x[j];
490 }
491 for (i=Mx; i<xs+xm+2; i++) {
492 for (j=0; j<dof; j++) x[i*dof+j] = x[(xs+xm-1)*dof+j];
493 }
494 }
495 for (i=xs-1; i<xs+xm+1; i++) {
496 struct _LimitInfo info;
497 PetscScalar *cjmpL,*cjmpR;
498 if ((i > sf-lsbwidth-2 && i < sf+1) || (i > fs-2 && i < fs+rsbwidth+1)) {
499 /* Determine the right eigenvectors R, where A = R \Lambda R^{-1} */
500 ierr = (*ctx->physics2.characteristic2)(ctx->physics2.user,dof,&x[i*dof],ctx->R,ctx->Rinv,ctx->speeds);CHKERRQ(ierr);
501 /* Evaluate jumps across interfaces (i-1, i) and (i, i+1), put in characteristic basis */
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.hxs = hxs;
517 info.hxf = hxf;
518 (*ctx->limit2)(&info,cjmpL,cjmpR,ctx->sf,ctx->fs,i,ctx->cslope);
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 uL = &ctx->uLR[0];
531 uR = &ctx->uLR[dof];
532 if (i == sf-lsbwidth) {
533 for (j=0; j<dof; j++) {
534 uL[j] = x[(i-1)*dof+j]+slope[(i-1)*dof+j]*hxs/2;
535 uR[j] = x[(i-0)*dof+j]-slope[(i-0)*dof+j]*hxs/2;
536 }
537 ierr = (*ctx->physics2.riemann2)(ctx->physics2.user,dof,uL,uR,ctx->flux,&maxspeed);CHKERRQ(ierr);
538 if (i < xs+xm) {
539 for (j=0; j<dof; j++) f[islow*dof+j] += ctx->flux[j]/hxs;
540 islow++;
541 }
542 }
543 if (i > sf-lsbwidth && i < sf) {
544 for (j=0; j<dof; j++) {
545 uL[j] = x[(i-1)*dof+j]+slope[(i-1)*dof+j]*hxs/2;
546 uR[j] = x[(i-0)*dof+j]-slope[(i-0)*dof+j]*hxs/2;
547 }
548 ierr = (*ctx->physics2.riemann2)(ctx->physics2.user,dof,uL,uR,ctx->flux,&maxspeed);CHKERRQ(ierr);
549 if (i > xs) {
550 for (j=0; j<dof; j++) f[(islow-1)*dof+j] -= ctx->flux[j]/hxs;
551 }
552 if (i < xs+xm) {
553 for (j=0; j<dof; j++) f[islow*dof+j] += ctx->flux[j]/hxs;
554 islow++;
555 }
556 }
557 if (i == sf) { /* interface between the slow region and the fast region */
558 for (j=0; j<dof; j++) {
559 uL[j] = x[(i-1)*dof+j]+slope[(i-1)*dof+j]*hxs/2;
560 uR[j] = x[(i-0)*dof+j]-slope[(i-0)*dof+j]*hxf/2;
561 }
562 ierr = (*ctx->physics2.riemann2)(ctx->physics2.user,dof,uL,uR,ctx->flux,&maxspeed);CHKERRQ(ierr);
563 if (i > xs) {
564 for (j=0; j<dof; j++) f[(islow-1)*dof+j] -= ctx->flux[j]/hxs;
565 }
566 }
567 if (i == fs) { /* interface between the fast region and the slow region */
568 for (j=0; j<dof; j++) {
569 uL[j] = x[(i-1)*dof+j]+slope[(i-1)*dof+j]*hxf/2;
570 uR[j] = x[(i-0)*dof+j]-slope[(i-0)*dof+j]*hxs/2;
571 }
572 ierr = (*ctx->physics2.riemann2)(ctx->physics2.user,dof,uL,uR,ctx->flux,&maxspeed);CHKERRQ(ierr);
573 if (i < xs+xm) {
574 for (j=0; j<dof; j++) f[islow*dof+j] += ctx->flux[j]/hxs;
575 islow++;
576 }
577 }
578 if (i > fs && i < fs+rsbwidth) {
579 for (j=0; j<dof; j++) {
580 uL[j] = x[(i-1)*dof+j]+slope[(i-1)*dof+j]*hxs/2;
581 uR[j] = x[(i-0)*dof+j]-slope[(i-0)*dof+j]*hxs/2;
582 }
583 ierr = (*ctx->physics2.riemann2)(ctx->physics2.user,dof,uL,uR,ctx->flux,&maxspeed);CHKERRQ(ierr);
584 if (i > xs) {
585 for (j=0; j<dof; j++) f[(islow-1)*dof+j] -= ctx->flux[j]/hxs;
586 }
587 if (i < xs+xm) {
588 for (j=0; j<dof; j++) f[islow*dof+j] += ctx->flux[j]/hxs;
589 islow++;
590 }
591 }
592 if (i == fs+rsbwidth) {
593 for (j=0; j<dof; j++) {
594 uL[j] = x[(i-1)*dof+j]+slope[(i-1)*dof+j]*hxs/2;
595 uR[j] = x[(i-0)*dof+j]-slope[(i-0)*dof+j]*hxs/2;
596 }
597 ierr = (*ctx->physics2.riemann2)(ctx->physics2.user,dof,uL,uR,ctx->flux,&maxspeed);CHKERRQ(ierr);
598 if (i > xs) {
599 for (j=0; j<dof; j++) f[(islow-1)*dof+j] -= ctx->flux[j]/hxs;
600 }
601 }
602 }
603 ierr = DMDAVecRestoreArray(da,Xloc,&x);CHKERRQ(ierr);
604 ierr = VecRestoreArray(F,&f);CHKERRQ(ierr);
605 ierr = DMDARestoreArray(da,PETSC_TRUE,&slope);CHKERRQ(ierr);
606 ierr = DMRestoreLocalVector(da,&Xloc);CHKERRQ(ierr);
607 PetscFunctionReturn(0);
608 }
609
610 /* --------------------------------- Finite Volume Solver for fast parts ----------------------------------- */
FVRHSFunctionfast_2WaySplit(TS ts,PetscReal time,Vec X,Vec F,void * vctx)611 PetscErrorCode FVRHSFunctionfast_2WaySplit(TS ts,PetscReal time,Vec X,Vec F,void *vctx)
612 {
613 FVCtx *ctx = (FVCtx*)vctx;
614 PetscErrorCode ierr;
615 PetscInt i,j,k,Mx,dof,xs,xm,ifast = 0,sf = ctx->sf,fs = ctx->fs;
616 PetscReal hxs,hxf;
617 PetscScalar *x,*f,*slope;
618 Vec Xloc;
619 DM da;
620
621 PetscFunctionBeginUser;
622 ierr = TSGetDM(ts,&da);CHKERRQ(ierr);
623 ierr = DMGetLocalVector(da,&Xloc);CHKERRQ(ierr);
624 ierr = DMDAGetInfo(da,0, &Mx,0,0, 0,0,0, &dof,0,0,0,0,0);CHKERRQ(ierr);
625 hxs = (ctx->xmax-ctx->xmin)*3.0/8.0/ctx->sf;
626 hxf = (ctx->xmax-ctx->xmin)/4.0/(ctx->fs-ctx->sf);
627 ierr = DMGlobalToLocalBegin(da,X,INSERT_VALUES,Xloc);CHKERRQ(ierr);
628 ierr = DMGlobalToLocalEnd(da,X,INSERT_VALUES,Xloc);CHKERRQ(ierr);
629 ierr = VecZeroEntries(F);CHKERRQ(ierr);
630 ierr = DMDAVecGetArray(da,Xloc,&x);CHKERRQ(ierr);
631 ierr = VecGetArray(F,&f);CHKERRQ(ierr);
632 ierr = DMDAGetArray(da,PETSC_TRUE,&slope);CHKERRQ(ierr);
633 ierr = DMDAGetCorners(da,&xs,0,0,&xm,0,0);CHKERRQ(ierr);
634
635 if (ctx->bctype == FVBC_OUTFLOW) {
636 for (i=xs-2; i<0; i++) {
637 for (j=0; j<dof; j++) x[i*dof+j] = x[j];
638 }
639 for (i=Mx; i<xs+xm+2; i++) {
640 for (j=0; j<dof; j++) x[i*dof+j] = x[(xs+xm-1)*dof+j];
641 }
642 }
643 for (i=xs-1; i<xs+xm+1; i++) {
644 struct _LimitInfo info;
645 PetscScalar *cjmpL,*cjmpR;
646 if (i > sf-2 && i < fs+1) {
647 ierr = (*ctx->physics2.characteristic2)(ctx->physics2.user,dof,&x[i*dof],ctx->R,ctx->Rinv,ctx->speeds);CHKERRQ(ierr);
648 ierr = PetscArrayzero(ctx->cjmpLR,2*dof);CHKERRQ(ierr);
649 cjmpL = &ctx->cjmpLR[0];
650 cjmpR = &ctx->cjmpLR[dof];
651 for (j=0; j<dof; j++) {
652 PetscScalar jmpL,jmpR;
653 jmpL = x[(i+0)*dof+j]-x[(i-1)*dof+j];
654 jmpR = x[(i+1)*dof+j]-x[(i+0)*dof+j];
655 for (k=0; k<dof; k++) {
656 cjmpL[k] += ctx->Rinv[k+j*dof]*jmpL;
657 cjmpR[k] += ctx->Rinv[k+j*dof]*jmpR;
658 }
659 }
660 /* Apply limiter to the left and right characteristic jumps */
661 info.m = dof;
662 info.hxs = hxs;
663 info.hxf = hxf;
664 (*ctx->limit2)(&info,cjmpL,cjmpR,ctx->sf,ctx->fs,i,ctx->cslope);
665 for (j=0; j<dof; j++) {
666 PetscScalar tmp = 0;
667 for (k=0; k<dof; k++) tmp += ctx->R[j+k*dof]*ctx->cslope[k];
668 slope[i*dof+j] = tmp;
669 }
670 }
671 }
672
673 for (i=xs; i<xs+xm+1; i++) {
674 PetscReal maxspeed;
675 PetscScalar *uL,*uR;
676 uL = &ctx->uLR[0];
677 uR = &ctx->uLR[dof];
678 if (i == sf) { /* interface between the slow region and the fast region */
679 for (j=0; j<dof; j++) {
680 uL[j] = x[(i-1)*dof+j]+slope[(i-1)*dof+j]*hxs/2;
681 uR[j] = x[(i-0)*dof+j]-slope[(i-0)*dof+j]*hxf/2;
682 }
683 ierr = (*ctx->physics2.riemann2)(ctx->physics2.user,dof,uL,uR,ctx->flux,&maxspeed);CHKERRQ(ierr);
684 if (i < xs+xm) {
685 for (j=0; j<dof; j++) f[ifast*dof+j] += ctx->flux[j]/hxf;
686 ifast++;
687 }
688 }
689 if (i > sf && i < fs) { /* fast region */
690 for (j=0; j<dof; j++) {
691 uL[j] = x[(i-1)*dof+j]+slope[(i-1)*dof+j]*hxf/2;
692 uR[j] = x[(i-0)*dof+j]-slope[(i-0)*dof+j]*hxf/2;
693 }
694 ierr = (*ctx->physics2.riemann2)(ctx->physics2.user,dof,uL,uR,ctx->flux,&maxspeed);CHKERRQ(ierr);
695 if (i > xs) {
696 for (j=0; j<dof; j++) f[(ifast-1)*dof+j] -= ctx->flux[j]/hxf;
697 }
698 if (i < xs+xm) {
699 for (j=0; j<dof; j++) f[ifast*dof+j] += ctx->flux[j]/hxf;
700 ifast++;
701 }
702 }
703 if (i == fs) { /* interface between the fast region and the slow region */
704 for (j=0; j<dof; j++) {
705 uL[j] = x[(i-1)*dof+j]+slope[(i-1)*dof+j]*hxf/2;
706 uR[j] = x[(i-0)*dof+j]-slope[(i-0)*dof+j]*hxs/2;
707 }
708 ierr = (*ctx->physics2.riemann2)(ctx->physics2.user,dof,uL,uR,ctx->flux,&maxspeed);CHKERRQ(ierr);
709 if (i > xs) {
710 for (j=0; j<dof; j++) f[(ifast-1)*dof+j] -= ctx->flux[j]/hxf;
711 }
712 }
713 }
714 ierr = DMDAVecRestoreArray(da,Xloc,&x);CHKERRQ(ierr);
715 ierr = VecRestoreArray(F,&f);CHKERRQ(ierr);
716 ierr = DMDARestoreArray(da,PETSC_TRUE,&slope);CHKERRQ(ierr);
717 ierr = DMRestoreLocalVector(da,&Xloc);CHKERRQ(ierr);
718 PetscFunctionReturn(0);
719 }
720
main(int argc,char * argv[])721 int main(int argc,char *argv[])
722 {
723 char lname[256] = "mc",physname[256] = "advect",final_fname[256] = "solution.m";
724 PetscFunctionList limiters = 0,physics = 0;
725 MPI_Comm comm;
726 TS ts;
727 DM da;
728 Vec X,X0,R;
729 FVCtx ctx;
730 PetscInt i,k,dof,xs,xm,Mx,draw = 0,count_slow,count_fast,islow = 0,ifast =0,islowbuffer = 0,*index_slow,*index_fast,*index_slowbuffer;
731 PetscBool view_final = PETSC_FALSE;
732 PetscReal ptime;
733 PetscErrorCode ierr;
734
735 ierr = PetscInitialize(&argc,&argv,0,help);if (ierr) return ierr;
736 comm = PETSC_COMM_WORLD;
737 ierr = PetscMemzero(&ctx,sizeof(ctx));CHKERRQ(ierr);
738
739 /* Register limiters to be available on the command line */
740 ierr = PetscFunctionListAdd(&limiters,"upwind" ,Limit2_Upwind);CHKERRQ(ierr);
741 ierr = PetscFunctionListAdd(&limiters,"lax-wendroff" ,Limit2_LaxWendroff);CHKERRQ(ierr);
742 ierr = PetscFunctionListAdd(&limiters,"beam-warming" ,Limit2_BeamWarming);CHKERRQ(ierr);
743 ierr = PetscFunctionListAdd(&limiters,"fromm" ,Limit2_Fromm);CHKERRQ(ierr);
744 ierr = PetscFunctionListAdd(&limiters,"minmod" ,Limit2_Minmod);CHKERRQ(ierr);
745 ierr = PetscFunctionListAdd(&limiters,"superbee" ,Limit2_Superbee);CHKERRQ(ierr);
746 ierr = PetscFunctionListAdd(&limiters,"mc" ,Limit2_MC);CHKERRQ(ierr);
747 ierr = PetscFunctionListAdd(&limiters,"koren3" ,Limit2_Koren3);CHKERRQ(ierr);
748
749 /* Register physical models to be available on the command line */
750 ierr = PetscFunctionListAdd(&physics,"advect" ,PhysicsCreate_Advect);CHKERRQ(ierr);
751
752 ctx.comm = comm;
753 ctx.cfl = 0.9;
754 ctx.bctype = FVBC_PERIODIC;
755 ctx.xmin = -1.0;
756 ctx.xmax = 1.0;
757 ierr = PetscOptionsBegin(comm,NULL,"Finite Volume solver options","");CHKERRQ(ierr);
758 ierr = PetscOptionsReal("-xmin","X min","",ctx.xmin,&ctx.xmin,NULL);CHKERRQ(ierr);
759 ierr = PetscOptionsReal("-xmax","X max","",ctx.xmax,&ctx.xmax,NULL);CHKERRQ(ierr);
760 ierr = PetscOptionsFList("-limit","Name of flux imiter to use","",limiters,lname,lname,sizeof(lname),NULL);CHKERRQ(ierr);
761 ierr = PetscOptionsInt("-draw","Draw solution vector, bitwise OR of (1=initial,2=final,4=final error)","",draw,&draw,NULL);CHKERRQ(ierr);
762 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);
763 ierr = PetscOptionsInt("-initial","Initial condition (depends on the physics)","",ctx.initial,&ctx.initial,NULL);CHKERRQ(ierr);
764 ierr = PetscOptionsBool("-exact","Compare errors with exact solution","",ctx.exact,&ctx.exact,NULL);CHKERRQ(ierr);
765 ierr = PetscOptionsBool("-simulation","Compare errors with reference solution","",ctx.simulation,&ctx.simulation,NULL);CHKERRQ(ierr);
766 ierr = PetscOptionsReal("-cfl","CFL number to time step at","",ctx.cfl,&ctx.cfl,NULL);CHKERRQ(ierr);
767 ierr = PetscOptionsEnum("-bc_type","Boundary condition","",FVBCTypes,(PetscEnum)ctx.bctype,(PetscEnum*)&ctx.bctype,NULL);CHKERRQ(ierr);
768 ierr = PetscOptionsInt("-hratio","Spacing ratio","",ctx.hratio,&ctx.hratio,NULL);CHKERRQ(ierr);
769 ierr = PetscOptionsEnd();CHKERRQ(ierr);
770
771 /* Choose the limiter from the list of registered limiters */
772 ierr = PetscFunctionListFind(limiters,lname,&ctx.limit2);CHKERRQ(ierr);
773 if (!ctx.limit2) SETERRQ1(PETSC_COMM_SELF,1,"Limiter '%s' not found",lname);
774
775 /* Choose the physics from the list of registered models */
776 {
777 PetscErrorCode (*r)(FVCtx*);
778 ierr = PetscFunctionListFind(physics,physname,&r);CHKERRQ(ierr);
779 if (!r) SETERRQ1(PETSC_COMM_SELF,1,"Physics '%s' not found",physname);
780 /* Create the physics, will set the number of fields and their names */
781 ierr = (*r)(&ctx);CHKERRQ(ierr);
782 }
783
784 /* Create a DMDA to manage the parallel grid */
785 ierr = DMDACreate1d(comm,DM_BOUNDARY_PERIODIC,50,ctx.physics2.dof,2,NULL,&da);CHKERRQ(ierr);
786 ierr = DMSetFromOptions(da);CHKERRQ(ierr);
787 ierr = DMSetUp(da);CHKERRQ(ierr);
788 /* Inform the DMDA of the field names provided by the physics. */
789 /* The names will be shown in the title bars when run with -ts_monitor_draw_solution */
790 for (i=0; i<ctx.physics2.dof; i++) {
791 ierr = DMDASetFieldName(da,i,ctx.physics2.fieldname[i]);CHKERRQ(ierr);
792 }
793 ierr = DMDAGetInfo(da,0, &Mx,0,0, 0,0,0, &dof,0,0,0,0,0);CHKERRQ(ierr);
794 ierr = DMDAGetCorners(da,&xs,0,0,&xm,0,0);CHKERRQ(ierr);
795
796 /* Set coordinates of cell centers */
797 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);
798
799 /* Allocate work space for the Finite Volume solver (so it doesn't have to be reallocated on each function evaluation) */
800 ierr = PetscMalloc4(dof*dof,&ctx.R,dof*dof,&ctx.Rinv,2*dof,&ctx.cjmpLR,1*dof,&ctx.cslope);CHKERRQ(ierr);
801 ierr = PetscMalloc3(2*dof,&ctx.uLR,dof,&ctx.flux,dof,&ctx.speeds);CHKERRQ(ierr);
802
803 /* Create a vector to store the solution and to save the initial state */
804 ierr = DMCreateGlobalVector(da,&X);CHKERRQ(ierr);
805 ierr = VecDuplicate(X,&X0);CHKERRQ(ierr);
806 ierr = VecDuplicate(X,&R);CHKERRQ(ierr);
807
808 /* create index for slow parts and fast parts,
809 count_slow + count_fast = Mx, counts_slow*hs = 0.5, counts_fast*hf = 0.5 */
810 count_slow = Mx/(1.0+ctx.hratio/3.0);
811 if (count_slow%2) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_USER,"Please adjust grid size Mx (-da_grid_x) and hratio (-hratio) so that Mx/(1+hartio/3) is even");
812 count_fast = Mx-count_slow;
813 ctx.sf = count_slow/2;
814 ctx.fs = ctx.sf+count_fast;
815 ierr = PetscMalloc1(xm*dof,&index_slow);CHKERRQ(ierr);
816 ierr = PetscMalloc1(xm*dof,&index_fast);CHKERRQ(ierr);
817 ierr = PetscMalloc1(6*dof,&index_slowbuffer);CHKERRQ(ierr);
818 if (((AdvectCtx*)ctx.physics2.user)->a > 0) {
819 ctx.lsbwidth = 2;
820 ctx.rsbwidth = 4;
821 } else {
822 ctx.lsbwidth = 4;
823 ctx.rsbwidth = 2;
824 }
825 for (i=xs; i<xs+xm; i++) {
826 if (i < ctx.sf-ctx.lsbwidth || i > ctx.fs+ctx.rsbwidth-1)
827 for (k=0; k<dof; k++) index_slow[islow++] = i*dof+k;
828 else if ((i >= ctx.sf-ctx.lsbwidth && i < ctx.sf) || (i > ctx.fs-1 && i <= ctx.fs+ctx.rsbwidth-1))
829 for (k=0; k<dof; k++) index_slowbuffer[islowbuffer++] = i*dof+k;
830 else
831 for (k=0; k<dof; k++) index_fast[ifast++] = i*dof+k;
832 }
833 ierr = ISCreateGeneral(PETSC_COMM_WORLD,islow,index_slow,PETSC_COPY_VALUES,&ctx.iss);CHKERRQ(ierr);
834 ierr = ISCreateGeneral(PETSC_COMM_WORLD,ifast,index_fast,PETSC_COPY_VALUES,&ctx.isf);CHKERRQ(ierr);
835 ierr = ISCreateGeneral(PETSC_COMM_WORLD,islowbuffer,index_slowbuffer,PETSC_COPY_VALUES,&ctx.issb);CHKERRQ(ierr);
836
837 /* Create a time-stepping object */
838 ierr = TSCreate(comm,&ts);CHKERRQ(ierr);
839 ierr = TSSetDM(ts,da);CHKERRQ(ierr);
840 ierr = TSSetRHSFunction(ts,R,FVRHSFunction_2WaySplit,&ctx);CHKERRQ(ierr);
841 ierr = TSRHSSplitSetIS(ts,"slow",ctx.iss);CHKERRQ(ierr);
842 ierr = TSRHSSplitSetIS(ts,"slowbuffer",ctx.issb);CHKERRQ(ierr);
843 ierr = TSRHSSplitSetIS(ts,"fast",ctx.isf);CHKERRQ(ierr);
844 ierr = TSRHSSplitSetRHSFunction(ts,"slow",NULL,FVRHSFunctionslow_2WaySplit,&ctx);CHKERRQ(ierr);
845 ierr = TSRHSSplitSetRHSFunction(ts,"fast",NULL,FVRHSFunctionfast_2WaySplit,&ctx);CHKERRQ(ierr);
846 ierr = TSRHSSplitSetRHSFunction(ts,"slowbuffer",NULL,FVRHSFunctionslowbuffer_2WaySplit,&ctx);CHKERRQ(ierr);
847
848 ierr = TSSetType(ts,TSSSP);CHKERRQ(ierr);
849 /*ierr = TSSetType(ts,TSMPRK);CHKERRQ(ierr);*/
850 ierr = TSSetMaxTime(ts,10);CHKERRQ(ierr);
851 ierr = TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);CHKERRQ(ierr);
852
853 /* Compute initial conditions and starting time step */
854 ierr = FVSample_2WaySplit(&ctx,da,0,X0);CHKERRQ(ierr);
855 ierr = FVRHSFunction_2WaySplit(ts,0,X0,X,(void*)&ctx);CHKERRQ(ierr); /* Initial function evaluation, only used to determine max speed */
856 ierr = VecCopy(X0,X);CHKERRQ(ierr); /* The function value was not used so we set X=X0 again */
857 ierr = TSSetTimeStep(ts,ctx.cfl/ctx.cfl_idt);CHKERRQ(ierr);
858 ierr = TSSetFromOptions(ts);CHKERRQ(ierr); /* Take runtime options */
859 ierr = SolutionStatsView(da,X,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
860 {
861 PetscInt steps;
862 PetscScalar mass_initial,mass_final,mass_difference,mass_differenceg;
863 const PetscScalar *ptr_X,*ptr_X0;
864 const PetscReal hs = (ctx.xmax-ctx.xmin)*3.0/4.0/count_slow;
865 const PetscReal hf = (ctx.xmax-ctx.xmin)/4.0/count_fast;
866
867 ierr = TSSolve(ts,X);CHKERRQ(ierr);
868 ierr = TSGetSolveTime(ts,&ptime);CHKERRQ(ierr);
869 ierr = TSGetStepNumber(ts,&steps);CHKERRQ(ierr);
870 /* calculate the total mass at initial time and final time */
871 mass_initial = 0.0;
872 mass_final = 0.0;
873 ierr = DMDAVecGetArrayRead(da,X0,(void*)&ptr_X0);CHKERRQ(ierr);
874 ierr = DMDAVecGetArrayRead(da,X,(void*)&ptr_X);CHKERRQ(ierr);
875 for (i=xs;i<xs+xm;i++) {
876 if (i < ctx.sf || i > ctx.fs-1) {
877 for (k=0; k<dof; k++) {
878 mass_initial = mass_initial + hs*ptr_X0[i*dof+k];
879 mass_final = mass_final + hs*ptr_X[i*dof+k];
880 }
881 } else {
882 for (k=0; k<dof; k++) {
883 mass_initial = mass_initial + hf*ptr_X0[i*dof+k];
884 mass_final = mass_final + hf*ptr_X[i*dof+k];
885 }
886 }
887 }
888 ierr = DMDAVecRestoreArrayRead(da,X0,(void*)&ptr_X0);CHKERRQ(ierr);
889 ierr = DMDAVecRestoreArrayRead(da,X,(void*)&ptr_X);CHKERRQ(ierr);
890 mass_difference = mass_final - mass_initial;
891 ierr = MPI_Allreduce(&mass_difference,&mass_differenceg,1,MPIU_SCALAR,MPIU_SUM,comm);CHKERRQ(ierr);
892 ierr = PetscPrintf(comm,"Mass difference %g\n",(double)mass_differenceg);CHKERRQ(ierr);
893 ierr = PetscPrintf(comm,"Final time %g, steps %D\n",(double)ptime,steps);CHKERRQ(ierr);
894 ierr = PetscPrintf(comm,"Maximum allowable stepsize according to CFL %g\n",(double)(1.0/ctx.cfl_idt));CHKERRQ(ierr);
895 if (ctx.exact) {
896 PetscReal nrm1=0;
897 ierr = SolutionErrorNorms_2WaySplit(&ctx,da,ptime,X,&nrm1);CHKERRQ(ierr);
898 ierr = PetscPrintf(comm,"Error ||x-x_e||_1 %g\n",(double)nrm1);CHKERRQ(ierr);
899 }
900 if (ctx.simulation) {
901 PetscReal nrm1=0;
902 PetscViewer fd;
903 char filename[PETSC_MAX_PATH_LEN] = "binaryoutput";
904 Vec XR;
905 PetscBool flg;
906 const PetscScalar *ptr_XR;
907 ierr = PetscOptionsGetString(NULL,NULL,"-f",filename,sizeof(filename),&flg);CHKERRQ(ierr);
908 if (!flg) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_USER,"Must indicate binary file with the -f option");
909 ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_READ,&fd);CHKERRQ(ierr);
910 ierr = VecDuplicate(X0,&XR);CHKERRQ(ierr);
911 ierr = VecLoad(XR,fd);CHKERRQ(ierr);
912 ierr = PetscViewerDestroy(&fd);CHKERRQ(ierr);
913 ierr = VecGetArrayRead(X,&ptr_X);CHKERRQ(ierr);
914 ierr = VecGetArrayRead(XR,&ptr_XR);CHKERRQ(ierr);
915 for (i=xs;i<xs+xm;i++) {
916 if (i < ctx.sf || i > ctx.fs-1)
917 for (k=0; k<dof; k++) nrm1 = nrm1 + hs*PetscAbs(ptr_X[i*dof+k]-ptr_XR[i*dof+k]);
918 else
919 for (k=0; k<dof; k++) nrm1 = nrm1 + hf*PetscAbs(ptr_X[i*dof+k]-ptr_XR[i*dof+k]);
920 }
921 ierr = VecRestoreArrayRead(X,&ptr_X);CHKERRQ(ierr);
922 ierr = VecRestoreArrayRead(XR,&ptr_XR);CHKERRQ(ierr);
923 ierr = PetscPrintf(comm,"Error ||x-x_e||_1 %g\n",(double)nrm1);CHKERRQ(ierr);
924 ierr = VecDestroy(&XR);CHKERRQ(ierr);
925 }
926 }
927
928 ierr = SolutionStatsView(da,X,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
929 if (draw & 0x1) {ierr = VecView(X0,PETSC_VIEWER_DRAW_WORLD);CHKERRQ(ierr);}
930 if (draw & 0x2) {ierr = VecView(X,PETSC_VIEWER_DRAW_WORLD);CHKERRQ(ierr);}
931 if (draw & 0x4) {
932 Vec Y;
933 ierr = VecDuplicate(X,&Y);CHKERRQ(ierr);
934 ierr = FVSample_2WaySplit(&ctx,da,ptime,Y);CHKERRQ(ierr);
935 ierr = VecAYPX(Y,-1,X);CHKERRQ(ierr);
936 ierr = VecView(Y,PETSC_VIEWER_DRAW_WORLD);CHKERRQ(ierr);
937 ierr = VecDestroy(&Y);CHKERRQ(ierr);
938 }
939
940 if (view_final) {
941 PetscViewer viewer;
942 ierr = PetscViewerASCIIOpen(PETSC_COMM_WORLD,final_fname,&viewer);CHKERRQ(ierr);
943 ierr = PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_MATLAB);CHKERRQ(ierr);
944 ierr = VecView(X,viewer);CHKERRQ(ierr);
945 ierr = PetscViewerPopFormat(viewer);CHKERRQ(ierr);
946 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);
947 }
948
949 /* Clean up */
950 ierr = (*ctx.physics2.destroy)(ctx.physics2.user);CHKERRQ(ierr);
951 for (i=0; i<ctx.physics2.dof; i++) {ierr = PetscFree(ctx.physics2.fieldname[i]);CHKERRQ(ierr);}
952 ierr = PetscFree4(ctx.R,ctx.Rinv,ctx.cjmpLR,ctx.cslope);CHKERRQ(ierr);
953 ierr = PetscFree3(ctx.uLR,ctx.flux,ctx.speeds);CHKERRQ(ierr);
954 ierr = VecDestroy(&X);CHKERRQ(ierr);
955 ierr = VecDestroy(&X0);CHKERRQ(ierr);
956 ierr = VecDestroy(&R);CHKERRQ(ierr);
957 ierr = DMDestroy(&da);CHKERRQ(ierr);
958 ierr = TSDestroy(&ts);CHKERRQ(ierr);
959 ierr = ISDestroy(&ctx.iss);CHKERRQ(ierr);
960 ierr = ISDestroy(&ctx.isf);CHKERRQ(ierr);
961 ierr = ISDestroy(&ctx.issb);CHKERRQ(ierr);
962 ierr = PetscFree(index_slow);CHKERRQ(ierr);
963 ierr = PetscFree(index_fast);CHKERRQ(ierr);
964 ierr = PetscFree(index_slowbuffer);CHKERRQ(ierr);
965 ierr = PetscFunctionListDestroy(&limiters);CHKERRQ(ierr);
966 ierr = PetscFunctionListDestroy(&physics);CHKERRQ(ierr);
967 ierr = PetscFinalize();
968 return ierr;
969 }
970
971 /*TEST
972
973 build:
974 requires: !complex
975 depends: finitevolume1d.c
976
977 test:
978 suffix: 1
979 args: -da_grid_x 60 -initial 7 -xmin -1 -xmax 1 -hratio 2 -limit mc -ts_dt 0.025 -ts_max_steps 24 -ts_type mprk -ts_mprk_type 2a22
980
981 test:
982 suffix: 2
983 args: -da_grid_x 60 -initial 7 -xmin -1 -xmax 1 -hratio 2 -limit mc -ts_dt 0.025 -ts_max_steps 24 -ts_type mprk -ts_mprk_type 2a22 -ts_use_splitrhsfunction 0
984 output_file: output/ex6_1.out
985
986 TEST*/
987