Actual source code: ex4.c

  1: /*
  2:   Note:
  3:     -hratio is the ratio between mesh size of coarse grids and fine grids
  4: */

  6: static const char help[] = "1D periodic Finite Volume solver in slope-limiter form with semidiscrete time stepping.\n"
  7:   "  advect      - Constant coefficient scalar advection\n"
  8:   "                u_t       + (a*u)_x               = 0\n"
  9:   "  shallow     - 1D Shallow water equations (Saint Venant System)\n"
 10:   "                h_t + (q)_x = 0 \n"
 11:   "                q_t + (\frac{q^2}{h} + g/2*h^2)_x = 0  \n"
 12:   "                where, h(x,t) denotes the height of the water and q(x,t) the momentum.\n"
 13:   "  for this toy problem, we choose different meshsizes for different sub-domains (slow-fast-slow), say\n"
 14:   "                hxs  = hratio*hxf \n"
 15:   "  where hxs and hxf are the grid spacings for coarse and fine grids respectively.\n"
 16:   "  exact       - Exact Riemann solver which usually needs to perform a Newton iteration to connect\n"
 17:   "                the states across shocks and rarefactions\n"
 18:   "  simulation  - use reference solution which is generated by smaller time step size to be true solution,\n"
 19:   "                also the reference solution should be generated by user and stored in a binary file.\n"
 20:   "  characteristic - Limit the characteristic variables, this is usually preferred (default)\n"
 21:   "  bc_type     - Boundary condition for the problem, options are: periodic, outflow, inflow "
 22:   "Several problem descriptions (initial data, physics specific features, boundary data) can be chosen with -initial N\n\n"
 23:   "The problem size should be set with -da_grid_x M\n\n";

 25: /*
 26:   Example:
 27:      ./ex4 -da_grid_x 40 -initial 1 -hratio 1 -limit mc -ts_dt 0.01 -ts_max_time 7.0 -ts_type mprk -ts_mprk_type 2a22 -ts_monitor_draw_solution -physics shallow -bc_type outflow -xmin 0 -xmax 50 -ts_use_splitrhsfunction 0
 28:      ./ex4 -da_grid_x 40 -initial 2 -hratio 1 -limit mc -ts_dt 0.01 -ts_max_time 2.5 -ts_type mprk -ts_mprk_type 2a22 -ts_monitor_draw_solution -physics shallow -bc_type outflow -xmin 0 -xmax 50 -ts_use_splitrhsfunction 0
 29:      ./ex4 -da_grid_x 40 -initial 3 -hratio 1 -limit mc -ts_dt 0.01 -ts_max_time 4.0 -ts_type mprk -ts_mprk_type 2a22 -ts_monitor_draw_solution -physics shallow -bc_type outflow -xmin 0 -xmax 50 -ts_use_splitrhsfunction 0
 30:      ./ex4 -da_grid_x 40 -initial 4 -hratio 1 -limit koren3 -ts_dt 0.01 -ts_max_time 4.0 -ts_type mprk -ts_mprk_type 2a22 -ts_monitor_draw_solution -physics shallow -bc_type outflow -xmin 0 -xmax 50 -ts_use_splitrhsfunction 0
 31:      ./ex4 -da_grid_x 40 -initial 5 -hratio 1 -limit mc -ts_dt 0.01 -ts_max_time 5.0 -ts_type mprk -ts_mprk_type 2a22 -ts_monitor_draw_solution -physics shallow -bc_type outflow -xmin 0 -xmax 50 -ts_use_splitrhsfunction 0

 33:   Contributed by: Aidan Hamilton <aidan@udel.edu>
 34: */

 36: #include <petscts.h>
 37: #include <petscdm.h>
 38: #include <petscdmda.h>
 39: #include <petscdraw.h>
 40: #include "finitevolume1d.h"
 41: #include <petsc/private/kernels/blockinvert.h>

 43: static inline PetscReal RangeMod(PetscReal a,PetscReal xmin,PetscReal xmax) { PetscReal range = xmax-xmin; return xmin +PetscFmodReal(range+PetscFmodReal(a,range),range); }
 44: static inline PetscReal MaxAbs(PetscReal a,PetscReal b) { return (PetscAbs(a) > PetscAbs(b)) ? a : b; }
 45: /* --------------------------------- Advection ----------------------------------- */
 46: typedef struct {
 47:   PetscReal a;                  /* advective velocity */
 48: } AdvectCtx;

 50: static PetscErrorCode PhysicsRiemann_Advect(void *vctx,PetscInt m,const PetscScalar *uL,const PetscScalar *uR,PetscScalar *flux,PetscReal *maxspeed)
 51: {
 52:   AdvectCtx *ctx = (AdvectCtx*)vctx;
 53:   PetscReal speed;

 56:   speed     = ctx->a;
 57:   flux[0]   = PetscMax(0,speed)*uL[0] + PetscMin(0,speed)*uR[0];
 58:   *maxspeed = speed;
 59:   return 0;
 60: }

 62: static PetscErrorCode PhysicsCharacteristic_Advect(void *vctx,PetscInt m,const PetscScalar *u,PetscScalar *X,PetscScalar *Xi,PetscReal *speeds)
 63: {
 64:   AdvectCtx *ctx = (AdvectCtx*)vctx;

 67:   X[0]      = 1.;
 68:   Xi[0]     = 1.;
 69:   speeds[0] = ctx->a;
 70:   return 0;
 71: }

 73: static PetscErrorCode PhysicsSample_Advect(void *vctx,PetscInt initial,FVBCType bctype,PetscReal xmin,PetscReal xmax,PetscReal t,PetscReal x,PetscReal *u)
 74: {
 75:   AdvectCtx *ctx = (AdvectCtx*)vctx;
 76:   PetscReal a    = ctx->a,x0;

 79:   switch (bctype) {
 80:     case FVBC_OUTFLOW:   x0 = x-a*t; break;
 81:     case FVBC_PERIODIC: x0 = RangeMod(x-a*t,xmin,xmax); break;
 82:     default: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"unknown BCType");
 83:   }
 84:   switch (initial) {
 85:     case 0: u[0] = (x0 < 0) ? 1 : -1; break;
 86:     case 1: u[0] = (x0 < 0) ? -1 : 1; break;
 87:     case 2: u[0] = (0 < x0 && x0 < 1) ? 1 : 0; break;
 88:     case 3: u[0] = PetscSinReal(2*PETSC_PI*x0); break;
 89:     case 4: u[0] = PetscAbs(x0); break;
 90:     case 5: u[0] = (x0 < 0 || x0 > 0.5) ? 0 : PetscSqr(PetscSinReal(2*PETSC_PI*x0)); break;
 91:     case 6: u[0] = (x0 < 0) ? 0 : ((x0 < 1) ? x0 : ((x0 < 2) ? 2-x0 : 0)); break;
 92:     case 7: u[0] = PetscPowReal(PetscSinReal(PETSC_PI*x0),10.0);break;
 93:     default: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"unknown initial condition");
 94:   }
 95:   return 0;
 96: }

 98: static PetscErrorCode PhysicsCreate_Advect(FVCtx *ctx)
 99: {
101:   AdvectCtx      *user;

104:   PetscNew(&user);
105:   ctx->physics2.sample2         = PhysicsSample_Advect;
106:   ctx->physics2.riemann2        = PhysicsRiemann_Advect;
107:   ctx->physics2.characteristic2 = PhysicsCharacteristic_Advect;
108:   ctx->physics2.destroy         = PhysicsDestroy_SimpleFree;
109:   ctx->physics2.user            = user;
110:   ctx->physics2.dof             = 1;

112:   PetscStrallocpy("u",&ctx->physics2.fieldname[0]);
113:   user->a = 1;
114:   PetscOptionsBegin(ctx->comm,ctx->prefix,"Options for advection","");
115:     PetscOptionsReal("-physics_advect_a","Speed","",user->a,&user->a,NULL);
116:   PetscOptionsEnd();
117:   return 0;
118: }
119: /* --------------------------------- Shallow Water ----------------------------------- */

121: typedef struct {
122:   PetscReal gravity;
123: } ShallowCtx;

125: static inline void ShallowFlux(ShallowCtx *phys,const PetscScalar *u,PetscScalar *f)
126: {
127:   f[0] = u[1];
128:   f[1] = PetscSqr(u[1])/u[0] + 0.5*phys->gravity*PetscSqr(u[0]);
129: }

131: static PetscErrorCode PhysicsRiemann_Shallow_Exact(void *vctx,PetscInt m,const PetscScalar *uL,const PetscScalar *uR,PetscScalar *flux,PetscReal *maxspeed)
132: {
133:   ShallowCtx                *phys = (ShallowCtx*)vctx;
134:   PetscScalar               g    = phys->gravity,ustar[2],cL,cR,c,cstar;
135:   struct {PetscScalar h,u;} L = {uL[0],uL[1]/uL[0]},R = {uR[0],uR[1]/uR[0]},star;
136:   PetscInt                  i;

140:   cL = PetscSqrtScalar(g*L.h);
141:   cR = PetscSqrtScalar(g*R.h);
142:   c  = PetscMax(cL,cR);
143:   {
144:     /* Solve for star state */
145:     const PetscInt maxits = 50;
146:     PetscScalar tmp,res,res0=0,h0,h = 0.5*(L.h + R.h); /* initial guess */
147:     h0 = h;
148:     for (i=0; i<maxits; i++) {
149:       PetscScalar fr,fl,dfr,dfl;
150:       fl = (L.h < h)
151:         ? PetscSqrtScalar(0.5*g*(h*h - L.h*L.h)*(1/L.h - 1/h)) /* shock */
152:         : 2*PetscSqrtScalar(g*h) - 2*PetscSqrtScalar(g*L.h);   /* rarefaction */
153:       fr = (R.h < h)
154:         ? PetscSqrtScalar(0.5*g*(h*h - R.h*R.h)*(1/R.h - 1/h)) /* shock */
155:         : 2*PetscSqrtScalar(g*h) - 2*PetscSqrtScalar(g*R.h);   /* rarefaction */
156:       res = R.u - L.u + fr + fl;
158:       if (PetscAbsScalar(res) < PETSC_SQRT_MACHINE_EPSILON || (i > 0 && PetscAbsScalar(h-h0) < PETSC_SQRT_MACHINE_EPSILON)) {
159:         star.h = h;
160:         star.u = L.u - fl;
161:         goto converged;
162:       } else if (i > 0 && PetscAbsScalar(res) >= PetscAbsScalar(res0)) {        /* Line search */
163:         h = 0.8*h0 + 0.2*h;
164:         continue;
165:       }
166:       /* Accept the last step and take another */
167:       res0 = res;
168:       h0 = h;
169:       dfl = (L.h < h) ? 0.5/fl*0.5*g*(-L.h*L.h/(h*h) - 1 + 2*h/L.h) : PetscSqrtScalar(g/h);
170:       dfr = (R.h < h) ? 0.5/fr*0.5*g*(-R.h*R.h/(h*h) - 1 + 2*h/R.h) : PetscSqrtScalar(g/h);
171:       tmp = h - res/(dfr+dfl);
172:       if (tmp <= 0) h /= 2;   /* Guard against Newton shooting off to a negative thickness */
173:       else h = tmp;
175:     }
176:     SETERRQ(PETSC_COMM_SELF,1,"Newton iteration for star.h diverged after %D iterations",i);
177:   }
178: converged:
179:   cstar = PetscSqrtScalar(g*star.h);
180:   if (L.u-cL < 0 && 0 < star.u-cstar) { /* 1-wave is sonic rarefaction */
181:     PetscScalar ufan[2];
182:     ufan[0] = 1/g*PetscSqr(L.u/3 + 2./3*cL);
183:     ufan[1] = PetscSqrtScalar(g*ufan[0])*ufan[0];
184:     ShallowFlux(phys,ufan,flux);
185:   } else if (star.u+cstar < 0 && 0 < R.u+cR) { /* 2-wave is sonic rarefaction */
186:     PetscScalar ufan[2];
187:     ufan[0] = 1/g*PetscSqr(R.u/3 - 2./3*cR);
188:     ufan[1] = -PetscSqrtScalar(g*ufan[0])*ufan[0];
189:     ShallowFlux(phys,ufan,flux);
190:   } else if ((L.h >= star.h && L.u-c >= 0) || (L.h<star.h && (star.h*star.u-L.h*L.u)/(star.h-L.h) > 0)) {
191:     /* 1-wave is right-travelling shock (supersonic) */
192:     ShallowFlux(phys,uL,flux);
193:   } else if ((star.h <= R.h && R.u+c <= 0) || (star.h>R.h && (R.h*R.u-star.h*star.h)/(R.h-star.h) < 0)) {
194:     /* 2-wave is left-travelling shock (supersonic) */
195:     ShallowFlux(phys,uR,flux);
196:   } else {
197:     ustar[0] = star.h;
198:     ustar[1] = star.h*star.u;
199:     ShallowFlux(phys,ustar,flux);
200:   }
201:   *maxspeed = MaxAbs(MaxAbs(star.u-cstar,star.u+cstar),MaxAbs(L.u-cL,R.u+cR));
202:   return 0;
203: }

205: static PetscErrorCode PhysicsRiemann_Shallow_Rusanov(void *vctx,PetscInt m,const PetscScalar *uL,const PetscScalar *uR,PetscScalar *flux,PetscReal *maxspeed)
206: {
207:   ShallowCtx                *phys = (ShallowCtx*)vctx;
208:   PetscScalar               g = phys->gravity,fL[2],fR[2],s;
209:   struct {PetscScalar h,u;} L = {uL[0],uL[1]/uL[0]},R = {uR[0],uR[1]/uR[0]};
210:   PetscReal                 tol = 1e-6;

213:   /* Positivity preserving modification*/
214:   if (L.h < tol) L.u = 0.0;
215:   if (R.h < tol) R.u = 0.0;

217:   /*simple pos preserve limiter*/
218:   if (L.h < 0) L.h = 0;
219:   if (R.h < 0) R.h = 0;

221:   ShallowFlux(phys,uL,fL);
222:   ShallowFlux(phys,uR,fR);

224:   s         = PetscMax(PetscAbs(L.u)+PetscSqrtScalar(g*L.h),PetscAbs(R.u)+PetscSqrtScalar(g*R.h));
225:   flux[0]   = 0.5*(fL[0] + fR[0]) + 0.5*s*(L.h - R.h);
226:   flux[1]   = 0.5*(fL[1] + fR[1]) + 0.5*s*(uL[1] - uR[1]);
227:   *maxspeed = s;
228:   return 0;
229: }

231: static PetscErrorCode PhysicsCharacteristic_Conservative(void *vctx,PetscInt m,const PetscScalar *u,PetscScalar *X,PetscScalar *Xi,PetscReal *speeds)
232: {
233:   PetscInt i,j;

236:   for (i=0; i<m; i++) {
237:     for (j=0; j<m; j++) Xi[i*m+j] = X[i*m+j] = (PetscScalar)(i==j);
238:     speeds[i] = PETSC_MAX_REAL; /* Indicates invalid */
239:   }
240:   return 0;
241: }

243: static PetscErrorCode PhysicsCharacteristic_Shallow(void *vctx,PetscInt m,const PetscScalar *u,PetscScalar *X,PetscScalar *Xi,PetscReal *speeds)
244: {
245:   ShallowCtx     *phys = (ShallowCtx*)vctx;
246:   PetscReal      c;
247:   PetscReal      tol = 1e-6;

250:   c           = PetscSqrtScalar(u[0]*phys->gravity);

252:   if (u[0] < tol) { /*Use conservative variables*/
253:     X[0*2+0]  = 1;
254:     X[0*2+1]  = 0;
255:     X[1*2+0]  = 0;
256:     X[1*2+1]  = 1;
257:   } else {
258:     speeds[0] = u[1]/u[0] - c;
259:     speeds[1] = u[1]/u[0] + c;
260:     X[0*2+0]  = 1;
261:     X[0*2+1]  = speeds[0];
262:     X[1*2+0]  = 1;
263:     X[1*2+1]  = speeds[1];
264:   }

266:   PetscArraycpy(Xi,X,4);
267:   PetscKernel_A_gets_inverse_A_2(Xi,0,PETSC_FALSE,NULL);
268:   return 0;
269: }

271: static PetscErrorCode PhysicsSample_Shallow(void *vctx,PetscInt initial,FVBCType bctype,PetscReal xmin,PetscReal xmax,PetscReal t,PetscReal x,PetscReal *u)
272: {
275:   switch (initial) {
276:     case 0:
277:       u[0] = (x < 0) ? 2 : 1; /* Standard Dam Break Problem */
278:       u[1] = (x < 0) ? 0 : 0;
279:       break;
280:     case 1:
281:       u[0] = (x < 10) ?   1 : 0.1; /*The Next 5 problems are standard Riemann problem tests */
282:       u[1] = (x < 10) ? 2.5 : 0;
283:       break;
284:     case 2:
285:       u[0] = (x < 25) ?  1 : 1;
286:       u[1] = (x < 25) ? -5 : 5;
287:       break;
288:     case 3:
289:       u[0] = (x < 20) ?  1 : 1e-12;
290:       u[1] = (x < 20) ?  0 : 0;
291:       break;
292:     case 4:
293:       u[0] = (x < 30) ? 1e-12 : 1;
294:       u[1] = (x < 30) ? 0 : 0;
295:       break;
296:     case 5:
297:       u[0] = (x < 25) ?  0.1 : 0.1;
298:       u[1] = (x < 25) ? -0.3 : 0.3;
299:       break;
300:     case 6:
301:       u[0] = 1+0.5*PetscSinReal(2*PETSC_PI*x);
302:       u[1] = 1*u[0];
303:       break;
304:     case 7:
305:       if (x < -0.1) {
306:        u[0] = 1e-9;
307:        u[1] = 0.0;
308:       } else if (x < 0.1) {
309:        u[0] = 1.0;
310:        u[1] = 0.0;
311:       } else {
312:        u[0] = 1e-9;
313:        u[1] = 0.0;
314:       }
315:       break;
316:     case 8:
317:      if (x < -0.1) {
318:        u[0] = 2;
319:        u[1] = 0.0;
320:       } else if (x < 0.1) {
321:        u[0] = 3.0;
322:        u[1] = 0.0;
323:       } else {
324:        u[0] = 2;
325:        u[1] = 0.0;
326:       }
327:       break;
328:     default: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"unknown initial condition");
329:   }
330:   return 0;
331: }

333: /* Implements inflow conditions for the given initial conditions. Which conditions are actually enforced depends on
334:    on the results of PhysicsSetInflowType_Shallow. */
335: static PetscErrorCode PhysicsInflow_Shallow(void *vctx,PetscReal t,PetscReal x,PetscReal *u)
336: {
337:   FVCtx          *ctx = (FVCtx*)vctx;

340:   if (ctx->bctype == FVBC_INFLOW) {
341:     switch (ctx->initial) {
342:       case 0:
343:       case 1:
344:       case 2:
345:       case 3:
346:       case 4:
347:       case 5:
348:         u[0] = 0; u[1] = 0.0; /* Left boundary conditions */
349:         u[2] = 0; u[3] = 0.0; /* Right boundary conditions */
350:         break;
351:       case 6:
352:         u[0] = 0; u[1] = 0.0; /* Left boundary conditions */
353:         u[2] = 0; u[3] = 0.0; /* Right boundary conditions */
354:         break;
355:       case 7:
356:         u[0] = 0; u[1] = 0.0; /* Left boundary conditions */
357:         u[2] = 0; u[3] = 0.0; /* Right boundary conditions */
358:         break;
359:       case 8:
360:         u[0] = 0; u[1] = 1.0; /* Left boundary conditions */
361:         u[2] = 0; u[3] = -1.0;/* Right boundary conditions */
362:         break;
363:       default: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"unknown initial condition");
364:     }
365:   }
366:   return 0;
367: }

369: /* Selects which boundary conditions are marked as inflow and which as outflow when FVBC_INFLOW is selected. */
370: static PetscErrorCode PhysicsSetInflowType_Shallow(FVCtx *ctx)
371: {
373:   switch (ctx->initial) {
374:     case 0:
375:     case 1:
376:     case 2:
377:     case 3:
378:     case 4:
379:     case 5:
380:     case 6:
381:     case 7:
382:     case 8: /* Fix left and right momentum, height is outflow */
383:       ctx->physics2.bcinflowindex[0] = PETSC_FALSE;
384:       ctx->physics2.bcinflowindex[1] = PETSC_TRUE;
385:       ctx->physics2.bcinflowindex[2] = PETSC_FALSE;
386:       ctx->physics2.bcinflowindex[3] = PETSC_TRUE;
387:       break;
388:     default: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"unknown initial condition");
389:   }
390:   return 0;
391: }

393: static PetscErrorCode PhysicsCreate_Shallow(FVCtx *ctx)
394: {
395:   PetscErrorCode    ierr;
396:   ShallowCtx        *user;
397:   PetscFunctionList rlist = 0,rclist = 0;
398:   char              rname[256] = "rusanov",rcname[256] = "characteristic";

401:   PetscNew(&user);
402:   ctx->physics2.sample2         = PhysicsSample_Shallow;
403:   ctx->physics2.inflow          = PhysicsInflow_Shallow;
404:   ctx->physics2.destroy         = PhysicsDestroy_SimpleFree;
405:   ctx->physics2.riemann2        = PhysicsRiemann_Shallow_Rusanov;
406:   ctx->physics2.characteristic2 = PhysicsCharacteristic_Shallow;
407:   ctx->physics2.user            = user;
408:   ctx->physics2.dof             = 2;

410:   PetscMalloc1(2*(ctx->physics2.dof),&ctx->physics2.bcinflowindex);
411:   PhysicsSetInflowType_Shallow(ctx);

413:   PetscStrallocpy("height",&ctx->physics2.fieldname[0]);
414:   PetscStrallocpy("momentum",&ctx->physics2.fieldname[1]);

416:   user->gravity = 9.81;

418:   RiemannListAdd_2WaySplit(&rlist,"exact",  PhysicsRiemann_Shallow_Exact);
419:   RiemannListAdd_2WaySplit(&rlist,"rusanov",PhysicsRiemann_Shallow_Rusanov);
420:   ReconstructListAdd_2WaySplit(&rclist,"characteristic",PhysicsCharacteristic_Shallow);
421:   ReconstructListAdd_2WaySplit(&rclist,"conservative",PhysicsCharacteristic_Conservative);
422:   PetscOptionsBegin(ctx->comm,ctx->prefix,"Options for Shallow","");
423:     PetscOptionsReal("-physics_shallow_gravity","Gravity","",user->gravity,&user->gravity,NULL);
424:     PetscOptionsFList("-physics_shallow_riemann","Riemann solver","",rlist,rname,rname,sizeof(rname),NULL);
425:     PetscOptionsFList("-physics_shallow_reconstruct","Reconstruction","",rclist,rcname,rcname,sizeof(rcname),NULL);
426:   PetscOptionsEnd();
427:   RiemannListFind_2WaySplit(rlist,rname,&ctx->physics2.riemann2);
428:   ReconstructListFind_2WaySplit(rclist,rcname,&ctx->physics2.characteristic2);
429:   PetscFunctionListDestroy(&rlist);
430:   PetscFunctionListDestroy(&rclist);
431:   return 0;
432: }

434: PetscErrorCode FVSample_2WaySplit(FVCtx *ctx,DM da,PetscReal time,Vec U)
435: {
436:   PetscScalar     *u,*uj,xj,xi;
437:   PetscInt        i,j,k,dof,xs,xm,Mx;
438:   const PetscInt  N = 200;
439:   PetscReal       hs,hf;

443:   DMDAGetInfo(da,0, &Mx,0,0, 0,0,0, &dof,0,0,0,0,0);
444:   DMDAGetCorners(da,&xs,0,0,&xm,0,0);
445:   DMDAVecGetArray(da,U,&u);
446:   PetscMalloc1(dof,&uj);
447:   hs   = (ctx->xmax-ctx->xmin)*3.0/8.0/ctx->sf;
448:   hf   = (ctx->xmax-ctx->xmin)/4.0/(ctx->fs-ctx->sf);
449:   for (i=xs; i<xs+xm; i++) {
450:     if (i < ctx->sf) {
451:       xi = ctx->xmin+0.5*hs+i*hs;
452:       /* Integrate over cell i using trapezoid rule with N points. */
453:       for (k=0; k<dof; k++) u[i*dof+k] = 0;
454:       for (j=0; j<N+1; j++) {
455:         xj = xi+hs*(j-N/2)/(PetscReal)N;
456:         (*ctx->physics2.sample2)(ctx->physics2.user,ctx->initial,ctx->bctype,ctx->xmin,ctx->xmax,time,xj,uj);
457:         for (k=0; k<dof; k++) u[i*dof+k] += ((j==0 || j==N) ? 0.5 : 1.0)*uj[k]/N;
458:       }
459:     } else if (i < ctx->fs) {
460:       xi = ctx->xmin+ctx->sf*hs+0.5*hf+(i-ctx->sf)*hf;
461:       /* Integrate over cell i using trapezoid rule with N points. */
462:       for (k=0; k<dof; k++) u[i*dof+k] = 0;
463:       for (j=0; j<N+1; j++) {
464:         xj = xi+hf*(j-N/2)/(PetscReal)N;
465:         (*ctx->physics2.sample2)(ctx->physics2.user,ctx->initial,ctx->bctype,ctx->xmin,ctx->xmax,time,xj,uj);
466:         for (k=0; k<dof; k++) u[i*dof+k] += ((j==0 || j==N) ? 0.5 : 1.0)*uj[k]/N;
467:       }
468:     } else {
469:       xi = ctx->xmin+ctx->sf*hs+(ctx->fs-ctx->sf)*hf+0.5*hs+(i-ctx->fs)*hs;
470:       /* Integrate over cell i using trapezoid rule with N points. */
471:       for (k=0; k<dof; k++) u[i*dof+k] = 0;
472:       for (j=0; j<N+1; j++) {
473:         xj = xi+hs*(j-N/2)/(PetscReal)N;
474:         (*ctx->physics2.sample2)(ctx->physics2.user,ctx->initial,ctx->bctype,ctx->xmin,ctx->xmax,time,xj,uj);
475:         for (k=0; k<dof; k++) u[i*dof+k] += ((j==0 || j==N) ? 0.5 : 1.0)*uj[k]/N;
476:       }
477:     }
478:   }
479:   DMDAVecRestoreArray(da,U,&u);
480:   PetscFree(uj);
481:   return 0;
482: }

484: static PetscErrorCode SolutionErrorNorms_2WaySplit(FVCtx *ctx,DM da,PetscReal t,Vec X,PetscReal *nrm1)
485: {
486:   Vec               Y;
487:   PetscInt          i,Mx;
488:   const PetscScalar *ptr_X,*ptr_Y;
489:   PetscReal         hs,hf;

492:   VecGetSize(X,&Mx);
493:   VecDuplicate(X,&Y);
494:   FVSample_2WaySplit(ctx,da,t,Y);
495:   hs   = (ctx->xmax-ctx->xmin)*3.0/8.0/ctx->sf;
496:   hf   = (ctx->xmax-ctx->xmin)/4.0/(ctx->fs-ctx->sf);
497:   VecGetArrayRead(X,&ptr_X);
498:   VecGetArrayRead(Y,&ptr_Y);
499:   for (i=0; i<Mx; i++) {
500:     if (i < ctx->sf || i > ctx->fs -1) *nrm1 +=  hs*PetscAbs(ptr_X[i]-ptr_Y[i]);
501:     else *nrm1 += hf*PetscAbs(ptr_X[i]-ptr_Y[i]);
502:   }
503:   VecRestoreArrayRead(X,&ptr_X);
504:   VecRestoreArrayRead(Y,&ptr_Y);
505:   VecDestroy(&Y);
506:   return 0;
507: }

509: PetscErrorCode FVRHSFunction_2WaySplit(TS ts,PetscReal time,Vec X,Vec F,void *vctx)
510: {
511:   FVCtx          *ctx = (FVCtx*)vctx;
512:   PetscInt       i,j,k,Mx,dof,xs,xm,sf = ctx->sf,fs = ctx->fs;
513:   PetscReal      hxf,hxs,cfl_idt = 0;
514:   PetscScalar    *x,*f,*slope;
515:   Vec            Xloc;
516:   DM             da;

519:   TSGetDM(ts,&da);
520:   DMGetLocalVector(da,&Xloc);                          /* Xloc contains ghost points                                     */
521:   DMDAGetInfo(da,0, &Mx,0,0, 0,0,0, &dof,0,0,0,0,0);   /* Mx is the number of center points                              */
522:   hxs  = (ctx->xmax-ctx->xmin)*3.0/8.0/ctx->sf;
523:   hxf  = (ctx->xmax-ctx->xmin)/4.0/(ctx->fs-ctx->sf);
524:   DMGlobalToLocalBegin(da,X,INSERT_VALUES,Xloc);       /* X is solution vector which does not contain ghost points       */
525:   DMGlobalToLocalEnd(da,X,INSERT_VALUES,Xloc);

527:   VecZeroEntries(F);                                   /* F is the right hand side function corresponds to center points */

529:   DMDAVecGetArray(da,Xloc,&x);
530:   DMDAVecGetArray(da,F,&f);
531:   DMDAGetArray(da,PETSC_TRUE,&slope);                  /* contains ghost points                                           */
532:   DMDAGetCorners(da,&xs,0,0,&xm,0,0);

534:   if (ctx->bctype == FVBC_OUTFLOW) {
535:     for (i=xs-2; i<0; i++) {
536:       for (j=0; j<dof; j++) x[i*dof+j] = x[j];
537:     }
538:     for (i=Mx; i<xs+xm+2; i++) {
539:       for (j=0; j<dof; j++) x[i*dof+j] = x[(xs+xm-1)*dof+j];
540:     }
541:   }

543:   if (ctx->bctype == FVBC_INFLOW) {
544:     /* See LeVeque, R. (2002). Finite Volume Methods for Hyperbolic Problems. doi:10.1017/CBO9780511791253
545:     pages 137-138 for the scheme. */
546:     if (xs==0) { /* Left Boundary */
547:       ctx->physics2.inflow(ctx,time,ctx->xmin,ctx->ub);
548:       for (j=0; j<dof; j++) {
549:         if (ctx->physics2.bcinflowindex[j]) {
550:           for (i=-2; i<0; i++) x[i*dof+j] = 2.0*ctx->ub[j]-x[-(i+1)*dof+j];
551:         } else {
552:           for (i=-2; i<0; i++) x[i*dof+j] = x[j]; /* Outflow */
553:         }
554:       }
555:     }
556:     if (xs+xm==Mx) { /* Right Boundary */
557:       ctx->physics2.inflow(ctx,time,ctx->xmax,ctx->ub);
558:       for (j=0; j<dof; j++) {
559:         if (ctx->physics2.bcinflowindex[dof+j]) {
560:           for (i=Mx; i<Mx+2; i++) x[i*dof+j] = 2.0*ctx->ub[dof+j]-x[(2*Mx-(i+1))*dof+j];
561:         } else {
562:           for (i=Mx; i<Mx+2; i++) x[i*dof+j] = x[(Mx-1)*dof+j]; /* Outflow */
563:         }
564:       }
565:     }
566:   }

568:   for (i=xs-1; i<xs+xm+1; i++) {
569:     struct _LimitInfo info;
570:     PetscScalar       *cjmpL,*cjmpR;
571:     /* Determine the right eigenvectors R, where A = R \Lambda R^{-1} */
572:     (*ctx->physics2.characteristic2)(ctx->physics2.user,dof,&x[i*dof],ctx->R,ctx->Rinv,ctx->speeds);
573:     /* Evaluate jumps across interfaces (i-1, i) and (i, i+1), put in characteristic basis */
574:     PetscArrayzero(ctx->cjmpLR,2*dof);
575:     cjmpL = &ctx->cjmpLR[0];
576:     cjmpR = &ctx->cjmpLR[dof];
577:     for (j=0; j<dof; j++) {
578:       PetscScalar jmpL,jmpR;
579:       jmpL = x[(i+0)*dof+j]-x[(i-1)*dof+j];
580:       jmpR = x[(i+1)*dof+j]-x[(i+0)*dof+j];
581:       for (k=0; k<dof; k++) {
582:         cjmpL[k] += ctx->Rinv[k+j*dof]*jmpL;
583:         cjmpR[k] += ctx->Rinv[k+j*dof]*jmpR;
584:       }
585:     }
586:     /* Apply limiter to the left and right characteristic jumps */
587:     info.m  = dof;
588:     info.hxs = hxs;
589:     info.hxf = hxf;
590:     (*ctx->limit2)(&info,cjmpL,cjmpR,ctx->sf,ctx->fs,i,ctx->cslope);
591:     for (j=0; j<dof; j++) {
592:       PetscScalar tmp = 0;
593:       for (k=0; k<dof; k++) tmp += ctx->R[j+k*dof]*ctx->cslope[k];
594:       slope[i*dof+j] = tmp;
595:     }
596:   }

598:   for (i=xs; i<xs+xm+1; i++) {
599:     PetscReal   maxspeed;
600:     PetscScalar *uL,*uR;
601:     uL = &ctx->uLR[0];
602:     uR = &ctx->uLR[dof];
603:     if (i < sf) { /* slow region */
604:       for (j=0; j<dof; j++) {
605:         uL[j] = x[(i-1)*dof+j]+slope[(i-1)*dof+j]*hxs/2;
606:         uR[j] = x[(i-0)*dof+j]-slope[(i-0)*dof+j]*hxs/2;
607:       }
608:       (*ctx->physics2.riemann2)(ctx->physics2.user,dof,uL,uR,ctx->flux,&maxspeed);
609:       if (i > xs) {
610:         for (j=0; j<dof; j++) f[(i-1)*dof+j] -= ctx->flux[j]/hxs;
611:       }
612:       if (i < xs+xm) {
613:         for (j=0; j<dof; j++) f[i*dof+j] += ctx->flux[j]/hxs;
614:       }
615:     } else if (i == sf) { /* interface between the slow region and the fast region */
616:       for (j=0; j<dof; j++) {
617:         uL[j] = x[(i-1)*dof+j]+slope[(i-1)*dof+j]*hxs/2;
618:         uR[j] = x[(i-0)*dof+j]-slope[(i-0)*dof+j]*hxf/2;
619:       }
620:       (*ctx->physics2.riemann2)(ctx->physics2.user,dof,uL,uR,ctx->flux,&maxspeed);
621:       if (i > xs) {
622:         for (j=0; j<dof; j++) f[(i-1)*dof+j] -= ctx->flux[j]/hxs;
623:       }
624:       if (i < xs+xm) {
625:         for (j=0; j<dof; j++) f[i*dof+j] += ctx->flux[j]/hxf;
626:       }
627:     } else if (i > sf && i < fs) { /* fast region */
628:       for (j=0; j<dof; j++) {
629:         uL[j] = x[(i-1)*dof+j]+slope[(i-1)*dof+j]*hxf/2;
630:         uR[j] = x[(i-0)*dof+j]-slope[(i-0)*dof+j]*hxf/2;
631:       }
632:       (*ctx->physics2.riemann2)(ctx->physics2.user,dof,uL,uR,ctx->flux,&maxspeed);
633:       if (i > xs) {
634:         for (j=0; j<dof; j++) f[(i-1)*dof+j] -= ctx->flux[j]/hxf;
635:       }
636:       if (i < xs+xm) {
637:         for (j=0; j<dof; j++) f[i*dof+j] += ctx->flux[j]/hxf;
638:       }
639:     } else if (i == fs) { /* interface between the fast region and the slow region */
640:       for (j=0; j<dof; j++) {
641:         uL[j] = x[(i-1)*dof+j]+slope[(i-1)*dof+j]*hxf/2;
642:         uR[j] = x[(i-0)*dof+j]-slope[(i-0)*dof+j]*hxs/2;
643:       }
644:       (*ctx->physics2.riemann2)(ctx->physics2.user,dof,uL,uR,ctx->flux,&maxspeed);
645:       if (i > xs) {
646:         for (j=0; j<dof; j++) f[(i-1)*dof+j] -= ctx->flux[j]/hxf;
647:       }
648:       if (i < xs+xm) {
649:         for (j=0; j<dof; j++) f[i*dof+j] += ctx->flux[j]/hxs;
650:       }
651:     } else { /* slow region */
652:       for (j=0; j<dof; j++) {
653:         uL[j] = x[(i-1)*dof+j]+slope[(i-1)*dof+j]*hxs/2;
654:         uR[j] = x[(i-0)*dof+j]-slope[(i-0)*dof+j]*hxs/2;
655:       }
656:       (*ctx->physics2.riemann2)(ctx->physics2.user,dof,uL,uR,ctx->flux,&maxspeed);
657:       cfl_idt = PetscMax(cfl_idt,PetscAbsScalar(maxspeed/hxs)); /* Max allowable value of 1/Delta t */
658:       if (i > xs) {
659:         for (j=0; j<dof; j++) f[(i-1)*dof+j] -= ctx->flux[j]/hxs;
660:       }
661:       if (i < xs+xm) {
662:         for (j=0; j<dof; j++) f[i*dof+j] += ctx->flux[j]/hxs;
663:       }
664:     }
665:   }
666:   DMDAVecRestoreArray(da,Xloc,&x);
667:   DMDAVecRestoreArray(da,F,&f);
668:   DMDARestoreArray(da,PETSC_TRUE,&slope);
669:   DMRestoreLocalVector(da,&Xloc);
670:   MPI_Allreduce(&cfl_idt,&ctx->cfl_idt,1,MPIU_SCALAR,MPIU_MAX,PetscObjectComm((PetscObject)da));
671:   if (0) {
672:     /* We need to a way to inform the TS of a CFL constraint, this is a debugging fragment */
673:     PetscReal dt,tnow;
674:     TSGetTimeStep(ts,&dt);
675:     TSGetTime(ts,&tnow);
676:     if (dt > 0.5/ctx->cfl_idt) {
677:       if (1) {
678:         PetscPrintf(ctx->comm,"Stability constraint exceeded at t=%g, dt %g > %g\n",(double)tnow,(double)dt,(double)(0.5/ctx->cfl_idt));
679:       } else SETERRQ(PETSC_COMM_SELF,1,"Stability constraint exceeded, %g > %g",(double)dt,(double)(ctx->cfl/ctx->cfl_idt));
680:     }
681:   }
682:   return 0;
683: }

685: /* --------------------------------- Finite Volume Solver for slow components ----------------------------------- */
686: PetscErrorCode FVRHSFunctionslow_2WaySplit(TS ts,PetscReal time,Vec X,Vec F,void *vctx)
687: {
688:   FVCtx          *ctx = (FVCtx*)vctx;
689:   PetscInt       i,j,k,Mx,dof,xs,xm,islow = 0,sf = ctx->sf,fs = ctx->fs,lsbwidth = ctx->lsbwidth,rsbwidth = ctx->rsbwidth;
690:   PetscReal      hxs,hxf,cfl_idt = 0;
691:   PetscScalar    *x,*f,*slope;
692:   Vec            Xloc;
693:   DM             da;

696:   TSGetDM(ts,&da);
697:   DMGetLocalVector(da,&Xloc);
698:   DMDAGetInfo(da,0, &Mx,0,0, 0,0,0, &dof,0,0,0,0,0);
699:   hxs  = (ctx->xmax-ctx->xmin)*3.0/8.0/ctx->sf;
700:   hxf  = (ctx->xmax-ctx->xmin)/4.0/(ctx->fs-ctx->sf);
701:   DMGlobalToLocalBegin(da,X,INSERT_VALUES,Xloc);
702:   DMGlobalToLocalEnd(da,X,INSERT_VALUES,Xloc);
703:   VecZeroEntries(F);
704:   DMDAVecGetArray(da,Xloc,&x);
705:   VecGetArray(F,&f);
706:   DMDAGetArray(da,PETSC_TRUE,&slope);
707:   DMDAGetCorners(da,&xs,0,0,&xm,0,0);

709:   if (ctx->bctype == FVBC_OUTFLOW) {
710:     for (i=xs-2; i<0; i++) {
711:       for (j=0; j<dof; j++) x[i*dof+j] = x[j];
712:     }
713:     for (i=Mx; i<xs+xm+2; i++) {
714:       for (j=0; j<dof; j++) x[i*dof+j] = x[(xs+xm-1)*dof+j];
715:     }
716:   }
717:   if (ctx->bctype == FVBC_INFLOW) {
718:     /* See LeVeque, R. (2002). Finite Volume Methods for Hyperbolic Problems. doi:10.1017/CBO9780511791253
719:     pages 137-138 for the scheme. */
720:     if (xs==0) { /* Left Boundary */
721:       ctx->physics2.inflow(ctx,time,ctx->xmin,ctx->ub);
722:       for (j=0; j<dof; j++) {
723:         if (ctx->physics2.bcinflowindex[j]==PETSC_TRUE) {
724:           for (i=-2; i<0; i++) x[i*dof+j] = 2.0*ctx->ub[j]-x[-(i+1)*dof+j];
725:         } else {
726:           for (i=-2; i<0; i++) x[i*dof+j] = x[j]; /* Outflow */
727:         }
728:       }
729:     }
730:     if (xs+xm==Mx) { /* Right Boundary */
731:       ctx->physics2.inflow(ctx,time,ctx->xmax,ctx->ub);
732:       for (j=0; j<dof; j++) {
733:         if (ctx->physics2.bcinflowindex[dof+j]==PETSC_TRUE) {
734:           for (i=Mx; i<Mx+2; i++) x[i*dof+j] = 2.0*ctx->ub[dof+j]-x[(2*Mx-(i+1))*dof+j];
735:         } else {
736:           for (i=Mx; i<Mx+2; i++) x[i*dof+j] = x[(Mx-1)*dof+j]; /* Outflow */
737:         }
738:       }
739:     }
740:   }
741:   for (i=xs-1; i<xs+xm+1; i++) {
742:     struct _LimitInfo info;
743:     PetscScalar       *cjmpL,*cjmpR;
744:     if (i < sf-lsbwidth+1 || i > fs+rsbwidth-2) { /* slow components and the first and last fast components */
745:       /* Determine the right eigenvectors R, where A = R \Lambda R^{-1} */
746:       (*ctx->physics2.characteristic2)(ctx->physics2.user,dof,&x[i*dof],ctx->R,ctx->Rinv,ctx->speeds);
747:       /* Evaluate jumps across interfaces (i-1, i) and (i, i+1), put in characteristic basis */
748:       PetscArrayzero(ctx->cjmpLR,2*dof);
749:       cjmpL = &ctx->cjmpLR[0];
750:       cjmpR = &ctx->cjmpLR[dof];
751:       for (j=0; j<dof; j++) {
752:         PetscScalar jmpL,jmpR;
753:         jmpL = x[(i+0)*dof+j]-x[(i-1)*dof+j];
754:         jmpR = x[(i+1)*dof+j]-x[(i+0)*dof+j];
755:         for (k=0; k<dof; k++) {
756:           cjmpL[k] += ctx->Rinv[k+j*dof]*jmpL;
757:           cjmpR[k] += ctx->Rinv[k+j*dof]*jmpR;
758:         }
759:       }
760:       /* Apply limiter to the left and right characteristic jumps */
761:       info.m  = dof;
762:       info.hxs = hxs;
763:       info.hxf = hxf;
764:       (*ctx->limit2)(&info,cjmpL,cjmpR,ctx->sf,ctx->fs,i,ctx->cslope);
765:       for (j=0; j<dof; j++) {
766:         PetscScalar tmp = 0;
767:         for (k=0; k<dof; k++) tmp += ctx->R[j+k*dof]*ctx->cslope[k];
768:           slope[i*dof+j] = tmp;
769:       }
770:     }
771:   }

773:   for (i=xs; i<xs+xm+1; i++) {
774:     PetscReal   maxspeed;
775:     PetscScalar *uL,*uR;
776:     uL = &ctx->uLR[0];
777:     uR = &ctx->uLR[dof];
778:     if (i < sf-lsbwidth) { /* slow region */
779:       for (j=0; j<dof; j++) {
780:         uL[j] = x[(i-1)*dof+j]+slope[(i-1)*dof+j]*hxs/2;
781:         uR[j] = x[(i-0)*dof+j]-slope[(i-0)*dof+j]*hxs/2;
782:       }
783:       (*ctx->physics2.riemann2)(ctx->physics2.user,dof,uL,uR,ctx->flux,&maxspeed);
784:       cfl_idt = PetscMax(cfl_idt,PetscAbsScalar(maxspeed/hxs)); /* Max allowable value of 1/Delta t */
785:       if (i > xs) {
786:         for (j=0; j<dof; j++) f[(islow-1)*dof+j] -= ctx->flux[j]/hxs;
787:       }
788:       if (i < xs+xm) {
789:         for (j=0; j<dof; j++) f[islow*dof+j] += ctx->flux[j]/hxs;
790:         islow++;
791:       }
792:     }
793:     if (i == sf-lsbwidth) { /* interface between the slow region and the fast region */
794:       for (j=0; j<dof; j++) {
795:         uL[j] = x[(i-1)*dof+j]+slope[(i-1)*dof+j]*hxs/2;
796:         uR[j] = x[(i-0)*dof+j]-slope[(i-0)*dof+j]*hxs/2;
797:       }
798:       (*ctx->physics2.riemann2)(ctx->physics2.user,dof,uL,uR,ctx->flux,&maxspeed);
799:       if (i > xs) {
800:         for (j=0; j<dof; j++) f[(islow-1)*dof+j] -= ctx->flux[j]/hxs;
801:       }
802:     }
803:     if (i == fs+rsbwidth) { /* slow region */
804:       for (j=0; j<dof; j++) {
805:         uL[j] = x[(i-1)*dof+j]+slope[(i-1)*dof+j]*hxs/2;
806:         uR[j] = x[(i-0)*dof+j]-slope[(i-0)*dof+j]*hxs/2;
807:       }
808:       (*ctx->physics2.riemann2)(ctx->physics2.user,dof,uL,uR,ctx->flux,&maxspeed);
809:       if (i < xs+xm) {
810:         for (j=0; j<dof; j++) f[islow*dof+j] += ctx->flux[j]/hxs;
811:         islow++;
812:       }
813:     }
814:     if (i > fs+rsbwidth) { /* slow region */
815:       for (j=0; j<dof; j++) {
816:         uL[j] = x[(i-1)*dof+j]+slope[(i-1)*dof+j]*hxs/2;
817:         uR[j] = x[(i-0)*dof+j]-slope[(i-0)*dof+j]*hxs/2;
818:       }
819:       (*ctx->physics2.riemann2)(ctx->physics2.user,dof,uL,uR,ctx->flux,&maxspeed);
820:       if (i > xs) {
821:         for (j=0; j<dof; j++) f[(islow-1)*dof+j] -= ctx->flux[j]/hxs;
822:       }
823:       if (i < xs+xm) {
824:         for (j=0; j<dof; j++) f[islow*dof+j] += ctx->flux[j]/hxs;
825:         islow++;
826:       }
827:     }
828:   }
829:   DMDAVecRestoreArray(da,Xloc,&x);
830:   VecRestoreArray(F,&f);
831:   DMDARestoreArray(da,PETSC_TRUE,&slope);
832:   DMRestoreLocalVector(da,&Xloc);
833:   MPI_Allreduce(&cfl_idt,&ctx->cfl_idt,1,MPIU_SCALAR,MPIU_MAX,PetscObjectComm((PetscObject)da));
834:   return 0;
835: }

837: PetscErrorCode FVRHSFunctionslowbuffer_2WaySplit(TS ts,PetscReal time,Vec X,Vec F,void *vctx)
838: {
839:   FVCtx          *ctx = (FVCtx*)vctx;
840:   PetscInt       i,j,k,Mx,dof,xs,xm,islow = 0,sf = ctx->sf,fs = ctx->fs,lsbwidth = ctx->lsbwidth,rsbwidth = ctx->rsbwidth;
841:   PetscReal      hxs,hxf;
842:   PetscScalar    *x,*f,*slope;
843:   Vec            Xloc;
844:   DM             da;

847:   TSGetDM(ts,&da);
848:   DMGetLocalVector(da,&Xloc);
849:   DMDAGetInfo(da,0, &Mx,0,0, 0,0,0, &dof,0,0,0,0,0);
850:   hxs  = (ctx->xmax-ctx->xmin)*3.0/8.0/ctx->sf;
851:   hxf  = (ctx->xmax-ctx->xmin)/4.0/(ctx->fs-ctx->sf);
852:   DMGlobalToLocalBegin(da,X,INSERT_VALUES,Xloc);
853:   DMGlobalToLocalEnd(da,X,INSERT_VALUES,Xloc);
854:   VecZeroEntries(F);
855:   DMDAVecGetArray(da,Xloc,&x);
856:   VecGetArray(F,&f);
857:   DMDAGetArray(da,PETSC_TRUE,&slope);
858:   DMDAGetCorners(da,&xs,0,0,&xm,0,0);

860:   if (ctx->bctype == FVBC_OUTFLOW) {
861:     for (i=xs-2; i<0; i++) {
862:       for (j=0; j<dof; j++) x[i*dof+j] = x[j];
863:     }
864:     for (i=Mx; i<xs+xm+2; i++) {
865:       for (j=0; j<dof; j++) x[i*dof+j] = x[(xs+xm-1)*dof+j];
866:     }
867:   }
868:   if (ctx->bctype == FVBC_INFLOW) {
869:     /* See LeVeque, R. (2002). Finite Volume Methods for Hyperbolic Problems. doi:10.1017/CBO9780511791253
870:     pages 137-138 for the scheme. */
871:     if (xs==0) { /* Left Boundary */
872:       ctx->physics2.inflow(ctx,time,ctx->xmin,ctx->ub);
873:       for (j=0; j<dof; j++) {
874:         if (ctx->physics2.bcinflowindex[j]==PETSC_TRUE) {
875:           for (i=-2; i<0; i++) x[i*dof+j] = 2.0*ctx->ub[j]-x[-(i+1)*dof+j];
876:         } else {
877:           for (i=-2; i<0; i++) x[i*dof+j] = x[j]; /* Outflow */
878:         }
879:       }
880:     }
881:     if (xs+xm==Mx) { /* Right Boundary */
882:       ctx->physics2.inflow(ctx,time,ctx->xmax,ctx->ub);
883:       for (j=0; j<dof; j++) {
884:         if (ctx->physics2.bcinflowindex[dof+j]==PETSC_TRUE) {
885:           for (i=Mx; i<Mx+2; i++) x[i*dof+j] = 2.0*ctx->ub[dof+j]-x[(2*Mx-(i+1))*dof+j];
886:         } else {
887:           for (i=Mx; i<Mx+2; i++) x[i*dof+j] = x[(Mx-1)*dof+j]; /* Outflow */
888:         }
889:       }
890:     }
891:   }
892:   for (i=xs-1; i<xs+xm+1; i++) {
893:     struct _LimitInfo info;
894:     PetscScalar       *cjmpL,*cjmpR;
895:     if ((i > sf-lsbwidth-2 && i < sf+1) || (i > fs-2 && i < fs+rsbwidth+1)) {
896:       /* Determine the right eigenvectors R, where A = R \Lambda R^{-1} */
897:       (*ctx->physics2.characteristic2)(ctx->physics2.user,dof,&x[i*dof],ctx->R,ctx->Rinv,ctx->speeds);
898:       /* Evaluate jumps across interfaces (i-1, i) and (i, i+1), put in characteristic basis */
899:       PetscArrayzero(ctx->cjmpLR,2*dof);
900:       cjmpL = &ctx->cjmpLR[0];
901:       cjmpR = &ctx->cjmpLR[dof];
902:       for (j=0; j<dof; j++) {
903:         PetscScalar jmpL,jmpR;
904:         jmpL = x[(i+0)*dof+j]-x[(i-1)*dof+j];
905:         jmpR = x[(i+1)*dof+j]-x[(i+0)*dof+j];
906:         for (k=0; k<dof; k++) {
907:           cjmpL[k] += ctx->Rinv[k+j*dof]*jmpL;
908:           cjmpR[k] += ctx->Rinv[k+j*dof]*jmpR;
909:         }
910:       }
911:       /* Apply limiter to the left and right characteristic jumps */
912:       info.m  = dof;
913:       info.hxs = hxs;
914:       info.hxf = hxf;
915:       (*ctx->limit2)(&info,cjmpL,cjmpR,ctx->sf,ctx->fs,i,ctx->cslope);
916:       for (j=0; j<dof; j++) {
917:         PetscScalar tmp = 0;
918:         for (k=0; k<dof; k++) tmp += ctx->R[j+k*dof]*ctx->cslope[k];
919:           slope[i*dof+j] = tmp;
920:       }
921:     }
922:   }

924:   for (i=xs; i<xs+xm+1; i++) {
925:     PetscReal   maxspeed;
926:     PetscScalar *uL,*uR;
927:     uL = &ctx->uLR[0];
928:     uR = &ctx->uLR[dof];
929:     if (i == sf-lsbwidth) {
930:       for (j=0; j<dof; j++) {
931:         uL[j] = x[(i-1)*dof+j]+slope[(i-1)*dof+j]*hxs/2;
932:         uR[j] = x[(i-0)*dof+j]-slope[(i-0)*dof+j]*hxs/2;
933:       }
934:       (*ctx->physics2.riemann2)(ctx->physics2.user,dof,uL,uR,ctx->flux,&maxspeed);
935:       if (i < xs+xm) {
936:         for (j=0; j<dof; j++) f[islow*dof+j] += ctx->flux[j]/hxs;
937:         islow++;
938:       }
939:     }
940:     if (i > sf-lsbwidth && i < sf) {
941:       for (j=0; j<dof; j++) {
942:         uL[j] = x[(i-1)*dof+j]+slope[(i-1)*dof+j]*hxs/2;
943:         uR[j] = x[(i-0)*dof+j]-slope[(i-0)*dof+j]*hxs/2;
944:       }
945:       (*ctx->physics2.riemann2)(ctx->physics2.user,dof,uL,uR,ctx->flux,&maxspeed);
946:       if (i > xs) {
947:         for (j=0; j<dof; j++) f[(islow-1)*dof+j] -= ctx->flux[j]/hxs;
948:       }
949:       if (i < xs+xm) {
950:         for (j=0; j<dof; j++) f[islow*dof+j] += ctx->flux[j]/hxs;
951:         islow++;
952:       }
953:     }
954:     if (i == sf) { /* interface between the slow region and the fast region */
955:       for (j=0; j<dof; j++) {
956:         uL[j] = x[(i-1)*dof+j]+slope[(i-1)*dof+j]*hxs/2;
957:         uR[j] = x[(i-0)*dof+j]-slope[(i-0)*dof+j]*hxf/2;
958:       }
959:       (*ctx->physics2.riemann2)(ctx->physics2.user,dof,uL,uR,ctx->flux,&maxspeed);
960:       if (i > xs) {
961:         for (j=0; j<dof; j++) f[(islow-1)*dof+j] -= ctx->flux[j]/hxs;
962:       }
963:     }
964:     if (i == fs) { /* interface between the fast region and the slow region */
965:       for (j=0; j<dof; j++) {
966:         uL[j] = x[(i-1)*dof+j]+slope[(i-1)*dof+j]*hxf/2;
967:         uR[j] = x[(i-0)*dof+j]-slope[(i-0)*dof+j]*hxs/2;
968:       }
969:       (*ctx->physics2.riemann2)(ctx->physics2.user,dof,uL,uR,ctx->flux,&maxspeed);
970:       if (i < xs+xm) {
971:         for (j=0; j<dof; j++) f[islow*dof+j] += ctx->flux[j]/hxs;
972:         islow++;
973:       }
974:     }
975:     if (i > fs && i < fs+rsbwidth) {
976:       for (j=0; j<dof; j++) {
977:         uL[j] = x[(i-1)*dof+j]+slope[(i-1)*dof+j]*hxs/2;
978:         uR[j] = x[(i-0)*dof+j]-slope[(i-0)*dof+j]*hxs/2;
979:       }
980:       (*ctx->physics2.riemann2)(ctx->physics2.user,dof,uL,uR,ctx->flux,&maxspeed);
981:       if (i > xs) {
982:         for (j=0; j<dof; j++) f[(islow-1)*dof+j] -= ctx->flux[j]/hxs;
983:       }
984:       if (i < xs+xm) {
985:         for (j=0; j<dof; j++) f[islow*dof+j] += ctx->flux[j]/hxs;
986:         islow++;
987:       }
988:     }
989:     if (i == fs+rsbwidth) {
990:       for (j=0; j<dof; j++) {
991:         uL[j] = x[(i-1)*dof+j]+slope[(i-1)*dof+j]*hxs/2;
992:         uR[j] = x[(i-0)*dof+j]-slope[(i-0)*dof+j]*hxs/2;
993:       }
994:       (*ctx->physics2.riemann2)(ctx->physics2.user,dof,uL,uR,ctx->flux,&maxspeed);
995:       if (i > xs) {
996:         for (j=0; j<dof; j++) f[(islow-1)*dof+j] -= ctx->flux[j]/hxs;
997:       }
998:     }
999:   }
1000:   DMDAVecRestoreArray(da,Xloc,&x);
1001:   VecRestoreArray(F,&f);
1002:   DMDARestoreArray(da,PETSC_TRUE,&slope);
1003:   DMRestoreLocalVector(da,&Xloc);
1004:   return 0;
1005: }

1007: /* --------------------------------- Finite Volume Solver for fast  parts ----------------------------------- */
1008: PetscErrorCode FVRHSFunctionfast_2WaySplit(TS ts,PetscReal time,Vec X,Vec F,void *vctx)
1009: {
1010:   FVCtx          *ctx = (FVCtx*)vctx;
1011:   PetscInt       i,j,k,Mx,dof,xs,xm,ifast = 0,sf = ctx->sf,fs = ctx->fs;
1012:   PetscReal      hxs,hxf;
1013:   PetscScalar    *x,*f,*slope;
1014:   Vec            Xloc;
1015:   DM             da;

1018:   TSGetDM(ts,&da);
1019:   DMGetLocalVector(da,&Xloc);
1020:   DMDAGetInfo(da,0, &Mx,0,0, 0,0,0, &dof,0,0,0,0,0);
1021:   hxs  = (ctx->xmax-ctx->xmin)*3.0/8.0/ctx->sf;
1022:   hxf  = (ctx->xmax-ctx->xmin)/4.0/(ctx->fs-ctx->sf);
1023:   DMGlobalToLocalBegin(da,X,INSERT_VALUES,Xloc);
1024:   DMGlobalToLocalEnd(da,X,INSERT_VALUES,Xloc);
1025:   VecZeroEntries(F);
1026:   DMDAVecGetArray(da,Xloc,&x);
1027:   VecGetArray(F,&f);
1028:   DMDAGetArray(da,PETSC_TRUE,&slope);
1029:   DMDAGetCorners(da,&xs,0,0,&xm,0,0);

1031:   if (ctx->bctype == FVBC_OUTFLOW) {
1032:     for (i=xs-2; i<0; i++) {
1033:       for (j=0; j<dof; j++) x[i*dof+j] = x[j];
1034:     }
1035:     for (i=Mx; i<xs+xm+2; i++) {
1036:       for (j=0; j<dof; j++) x[i*dof+j] = x[(xs+xm-1)*dof+j];
1037:     }
1038:   }
1039:   if (ctx->bctype == FVBC_INFLOW) {
1040:     /* See LeVeque, R. (2002). Finite Volume Methods for Hyperbolic Problems. doi:10.1017/CBO9780511791253
1041:     pages 137-138 for the scheme. */
1042:     if (xs==0) { /* Left Boundary */
1043:       ctx->physics2.inflow(ctx,time,ctx->xmin,ctx->ub);
1044:       for (j=0; j<dof; j++) {
1045:         if (ctx->physics2.bcinflowindex[j]==PETSC_TRUE) {
1046:           for (i=-2; i<0; i++) x[i*dof+j] = 2.0*ctx->ub[j]-x[-(i+1)*dof+j];
1047:         } else {
1048:           for (i=-2; i<0; i++) x[i*dof+j] = x[j]; /* Outflow */
1049:         }
1050:       }
1051:     }
1052:     if (xs+xm==Mx) { /* Right Boundary */
1053:       ctx->physics2.inflow(ctx,time,ctx->xmax,ctx->ub);
1054:       for (j=0; j<dof; j++) {
1055:         if (ctx->physics2.bcinflowindex[dof+j]==PETSC_TRUE) {
1056:           for (i=Mx; i<Mx+2; i++) x[i*dof+j] = 2.0*ctx->ub[dof+j]-x[(2*Mx-(i+1))*dof+j];
1057:         } else {
1058:           for (i=Mx; i<Mx+2; i++) x[i*dof+j] = x[(Mx-1)*dof+j]; /* Outflow */
1059:         }
1060:       }
1061:     }
1062:   }
1063:   for (i=xs-1; i<xs+xm+1; i++) {
1064:     struct _LimitInfo info;
1065:     PetscScalar       *cjmpL,*cjmpR;
1066:     if (i > sf-2 && i < fs+1) {
1067:       (*ctx->physics2.characteristic2)(ctx->physics2.user,dof,&x[i*dof],ctx->R,ctx->Rinv,ctx->speeds);
1068:       PetscArrayzero(ctx->cjmpLR,2*dof);
1069:       cjmpL = &ctx->cjmpLR[0];
1070:       cjmpR = &ctx->cjmpLR[dof];
1071:       for (j=0; j<dof; j++) {
1072:         PetscScalar jmpL,jmpR;
1073:         jmpL = x[(i+0)*dof+j]-x[(i-1)*dof+j];
1074:         jmpR = x[(i+1)*dof+j]-x[(i+0)*dof+j];
1075:         for (k=0; k<dof; k++) {
1076:           cjmpL[k] += ctx->Rinv[k+j*dof]*jmpL;
1077:           cjmpR[k] += ctx->Rinv[k+j*dof]*jmpR;
1078:         }
1079:       }
1080:       /* Apply limiter to the left and right characteristic jumps */
1081:       info.m  = dof;
1082:       info.hxs = hxs;
1083:       info.hxf = hxf;
1084:       (*ctx->limit2)(&info,cjmpL,cjmpR,ctx->sf,ctx->fs,i,ctx->cslope);
1085:       for (j=0; j<dof; j++) {
1086:       PetscScalar tmp = 0;
1087:       for (k=0; k<dof; k++) tmp += ctx->R[j+k*dof]*ctx->cslope[k];
1088:         slope[i*dof+j] = tmp;
1089:       }
1090:     }
1091:   }

1093:   for (i=xs; i<xs+xm+1; i++) {
1094:     PetscReal   maxspeed;
1095:     PetscScalar *uL,*uR;
1096:     uL = &ctx->uLR[0];
1097:     uR = &ctx->uLR[dof];
1098:     if (i == sf) { /* interface between the slow region and the fast region */
1099:       for (j=0; j<dof; j++) {
1100:         uL[j] = x[(i-1)*dof+j]+slope[(i-1)*dof+j]*hxs/2;
1101:         uR[j] = x[(i-0)*dof+j]-slope[(i-0)*dof+j]*hxf/2;
1102:       }
1103:       (*ctx->physics2.riemann2)(ctx->physics2.user,dof,uL,uR,ctx->flux,&maxspeed);
1104:       if (i < xs+xm) {
1105:         for (j=0; j<dof; j++) f[ifast*dof+j] += ctx->flux[j]/hxf;
1106:         ifast++;
1107:       }
1108:     }
1109:     if (i > sf && i < fs) { /* fast region */
1110:       for (j=0; j<dof; j++) {
1111:         uL[j] = x[(i-1)*dof+j]+slope[(i-1)*dof+j]*hxf/2;
1112:         uR[j] = x[(i-0)*dof+j]-slope[(i-0)*dof+j]*hxf/2;
1113:       }
1114:       (*ctx->physics2.riemann2)(ctx->physics2.user,dof,uL,uR,ctx->flux,&maxspeed);
1115:       if (i > xs) {
1116:         for (j=0; j<dof; j++) f[(ifast-1)*dof+j] -= ctx->flux[j]/hxf;
1117:       }
1118:       if (i < xs+xm) {
1119:         for (j=0; j<dof; j++) f[ifast*dof+j] += ctx->flux[j]/hxf;
1120:         ifast++;
1121:       }
1122:     }
1123:     if (i == fs) { /* interface between the fast region and the slow region */
1124:       for (j=0; j<dof; j++) {
1125:         uL[j] = x[(i-1)*dof+j]+slope[(i-1)*dof+j]*hxf/2;
1126:         uR[j] = x[(i-0)*dof+j]-slope[(i-0)*dof+j]*hxs/2;
1127:       }
1128:       (*ctx->physics2.riemann2)(ctx->physics2.user,dof,uL,uR,ctx->flux,&maxspeed);
1129:       if (i > xs) {
1130:         for (j=0; j<dof; j++) f[(ifast-1)*dof+j] -= ctx->flux[j]/hxf;
1131:       }
1132:     }
1133:   }
1134:   DMDAVecRestoreArray(da,Xloc,&x);
1135:   VecRestoreArray(F,&f);
1136:   DMDARestoreArray(da,PETSC_TRUE,&slope);
1137:   DMRestoreLocalVector(da,&Xloc);
1138:   return 0;
1139: }

1141: int main(int argc,char *argv[])
1142: {
1143:   char              lname[256] = "mc",physname[256] = "advect",final_fname[256] = "solution.m";
1144:   PetscFunctionList limiters   = 0,physics = 0;
1145:   MPI_Comm          comm;
1146:   TS                ts;
1147:   DM                da;
1148:   Vec               X,X0,R;
1149:   FVCtx             ctx;
1150:   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;
1151:   PetscBool         view_final = PETSC_FALSE;
1152:   PetscReal         ptime,maxtime;
1153:   PetscErrorCode    ierr;

1155:   PetscInitialize(&argc,&argv,0,help);
1156:   comm = PETSC_COMM_WORLD;
1157:   PetscMemzero(&ctx,sizeof(ctx));

1159:   /* Register limiters to be available on the command line */
1160:   PetscFunctionListAdd(&limiters,"upwind"              ,Limit2_Upwind);
1161:   PetscFunctionListAdd(&limiters,"lax-wendroff"        ,Limit2_LaxWendroff);
1162:   PetscFunctionListAdd(&limiters,"beam-warming"        ,Limit2_BeamWarming);
1163:   PetscFunctionListAdd(&limiters,"fromm"               ,Limit2_Fromm);
1164:   PetscFunctionListAdd(&limiters,"minmod"              ,Limit2_Minmod);
1165:   PetscFunctionListAdd(&limiters,"superbee"            ,Limit2_Superbee);
1166:   PetscFunctionListAdd(&limiters,"mc"                  ,Limit2_MC);
1167:   PetscFunctionListAdd(&limiters,"koren3"              ,Limit2_Koren3);

1169:   /* Register physical models to be available on the command line */
1170:   PetscFunctionListAdd(&physics,"shallow"         ,PhysicsCreate_Shallow);
1171:   PetscFunctionListAdd(&physics,"advect"          ,PhysicsCreate_Advect);

1173:   ctx.comm    = comm;
1174:   ctx.cfl     = 0.9;
1175:   ctx.bctype  = FVBC_PERIODIC;
1176:   ctx.xmin    = -1.0;
1177:   ctx.xmax    = 1.0;
1178:   ctx.initial = 1;
1179:   ctx.hratio  = 2;
1180:   maxtime     = 10.0;
1181:   ctx.simulation = PETSC_FALSE;
1182:   PetscOptionsBegin(comm,NULL,"Finite Volume solver options","");
1183:   PetscOptionsReal("-xmin","X min","",ctx.xmin,&ctx.xmin,NULL);
1184:   PetscOptionsReal("-xmax","X max","",ctx.xmax,&ctx.xmax,NULL);
1185:   PetscOptionsFList("-limit","Name of flux imiter to use","",limiters,lname,lname,sizeof(lname),NULL);
1186:   PetscOptionsFList("-physics","Name of physics model to use","",physics,physname,physname,sizeof(physname),NULL);
1187:   PetscOptionsInt("-draw","Draw solution vector, bitwise OR of (1=initial,2=final,4=final error)","",draw,&draw,NULL);
1188:   PetscOptionsString("-view_final","Write final solution in ASCII MATLAB format to given file name","",final_fname,final_fname,sizeof(final_fname),&view_final);
1189:   PetscOptionsInt("-initial","Initial condition (depends on the physics)","",ctx.initial,&ctx.initial,NULL);
1190:   PetscOptionsBool("-exact","Compare errors with exact solution","",ctx.exact,&ctx.exact,NULL);
1191:   PetscOptionsBool("-simulation","Compare errors with reference solution","",ctx.simulation,&ctx.simulation,NULL);
1192:   PetscOptionsReal("-cfl","CFL number to time step at","",ctx.cfl,&ctx.cfl,NULL);
1193:   PetscOptionsEnum("-bc_type","Boundary condition","",FVBCTypes,(PetscEnum)ctx.bctype,(PetscEnum*)&ctx.bctype,NULL);
1194:   PetscOptionsInt("-hratio","Spacing ratio","",ctx.hratio,&ctx.hratio,NULL);
1195:   PetscOptionsEnd();

1197:   /* Choose the limiter from the list of registered limiters */
1198:   PetscFunctionListFind(limiters,lname,&ctx.limit2);

1201:   /* Choose the physics from the list of registered models */
1202:   {
1203:     PetscErrorCode (*r)(FVCtx*);
1204:     PetscFunctionListFind(physics,physname,&r);
1206:     /* Create the physics, will set the number of fields and their names */
1207:     (*r)(&ctx);
1208:   }

1210:   /* Create a DMDA to manage the parallel grid */
1211:   DMDACreate1d(comm,DM_BOUNDARY_PERIODIC,50,ctx.physics2.dof,2,NULL,&da);
1212:   DMSetFromOptions(da);
1213:   DMSetUp(da);
1214:   /* Inform the DMDA of the field names provided by the physics. */
1215:   /* The names will be shown in the title bars when run with -ts_monitor_draw_solution */
1216:   for (i=0; i<ctx.physics2.dof; i++) {
1217:     DMDASetFieldName(da,i,ctx.physics2.fieldname[i]);
1218:   }
1219:   DMDAGetInfo(da,0, &Mx,0,0, 0,0,0, &dof,0,0,0,0,0);
1220:   DMDAGetCorners(da,&xs,0,0,&xm,0,0);

1222:   /* Set coordinates of cell centers */
1223:   DMDASetUniformCoordinates(da,ctx.xmin+0.5*(ctx.xmax-ctx.xmin)/Mx,ctx.xmax+0.5*(ctx.xmax-ctx.xmin)/Mx,0,0,0,0);

1225:   /* Allocate work space for the Finite Volume solver (so it doesn't have to be reallocated on each function evaluation) */
1226:   PetscMalloc4(dof*dof,&ctx.R,dof*dof,&ctx.Rinv,2*dof,&ctx.cjmpLR,1*dof,&ctx.cslope);
1227:   PetscMalloc3(2*dof,&ctx.uLR,dof,&ctx.flux,dof,&ctx.speeds);
1228:   PetscMalloc1(2*dof,&ctx.ub);

1230:   /* Create a vector to store the solution and to save the initial state */
1231:   DMCreateGlobalVector(da,&X);
1232:   VecDuplicate(X,&X0);
1233:   VecDuplicate(X,&R);

1235:   /* create index for slow parts and fast parts,
1236:      count_slow + count_fast = Mx, counts_slow*hs = 0.5, counts_fast*hf = 0.5 */
1237:   count_slow = Mx/(1.0+ctx.hratio/3.0);
1239:   count_fast = Mx-count_slow;
1240:   ctx.sf = count_slow/2;
1241:   ctx.fs = ctx.sf+count_fast;
1242:   PetscMalloc1(xm*dof,&index_slow);
1243:   PetscMalloc1(xm*dof,&index_fast);
1244:   PetscMalloc1(8*dof,&index_slowbuffer);
1245:   ctx.lsbwidth = 4;
1246:   ctx.rsbwidth = 4;

1248:   for (i=xs; i<xs+xm; i++) {
1249:     if (i < ctx.sf-ctx.lsbwidth || i > ctx.fs+ctx.rsbwidth-1)
1250:       for (k=0; k<dof; k++) index_slow[islow++] = i*dof+k;
1251:     else if ((i >= ctx.sf-ctx.lsbwidth && i < ctx.sf) || (i > ctx.fs-1 && i <= ctx.fs+ctx.rsbwidth-1))
1252:       for (k=0; k<dof; k++) index_slowbuffer[islowbuffer++] = i*dof+k;
1253:     else
1254:       for (k=0; k<dof; k++) index_fast[ifast++] = i*dof+k;
1255:   }
1256:   ISCreateGeneral(PETSC_COMM_WORLD,islow,index_slow,PETSC_COPY_VALUES,&ctx.iss);
1257:   ISCreateGeneral(PETSC_COMM_WORLD,ifast,index_fast,PETSC_COPY_VALUES,&ctx.isf);
1258:   ISCreateGeneral(PETSC_COMM_WORLD,islowbuffer,index_slowbuffer,PETSC_COPY_VALUES,&ctx.issb);

1260:   /* Create a time-stepping object */
1261:   TSCreate(comm,&ts);
1262:   TSSetDM(ts,da);
1263:   TSSetRHSFunction(ts,R,FVRHSFunction_2WaySplit,&ctx);
1264:   TSRHSSplitSetIS(ts,"slow",ctx.iss);
1265:   TSRHSSplitSetIS(ts,"slowbuffer",ctx.issb);
1266:   TSRHSSplitSetIS(ts,"fast",ctx.isf);
1267:   TSRHSSplitSetRHSFunction(ts,"slow",NULL,FVRHSFunctionslow_2WaySplit,&ctx);
1268:   TSRHSSplitSetRHSFunction(ts,"fast",NULL,FVRHSFunctionfast_2WaySplit,&ctx);
1269:   TSRHSSplitSetRHSFunction(ts,"slowbuffer",NULL,FVRHSFunctionslowbuffer_2WaySplit,&ctx);

1271:   TSSetType(ts,TSMPRK);
1272:   TSSetMaxTime(ts,maxtime);
1273:   TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);

1275:   /* Compute initial conditions and starting time step */
1276:   FVSample_2WaySplit(&ctx,da,0,X0);
1277:   FVRHSFunction_2WaySplit(ts,0,X0,X,(void*)&ctx); /* Initial function evaluation, only used to determine max speed */
1278:   VecCopy(X0,X);                        /* The function value was not used so we set X=X0 again */
1279:   TSSetTimeStep(ts,ctx.cfl/ctx.cfl_idt);
1280:   TSSetFromOptions(ts); /* Take runtime options */
1281:   SolutionStatsView(da,X,PETSC_VIEWER_STDOUT_WORLD);
1282:   {
1283:     PetscInt          steps;
1284:     PetscScalar       mass_initial,mass_final,mass_difference,mass_differenceg;
1285:     const PetscScalar *ptr_X,*ptr_X0;
1286:     const PetscReal   hs = (ctx.xmax-ctx.xmin)*3.0/4.0/count_slow;
1287:     const PetscReal   hf = (ctx.xmax-ctx.xmin)/4.0/count_fast;

1289:     TSSolve(ts,X);
1290:     TSGetSolveTime(ts,&ptime);
1291:     TSGetStepNumber(ts,&steps);
1292:     /* calculate the total mass at initial time and final time */
1293:     mass_initial = 0.0;
1294:     mass_final   = 0.0;
1295:     DMDAVecGetArrayRead(da,X0,(void*)&ptr_X0);
1296:     DMDAVecGetArrayRead(da,X,(void*)&ptr_X);
1297:     for (i=xs;i<xs+xm;i++) {
1298:       if (i < ctx.sf || i > ctx.fs-1) {
1299:         for (k=0; k<dof; k++) {
1300:           mass_initial = mass_initial + hs*ptr_X0[i*dof+k];
1301:           mass_final = mass_final + hs*ptr_X[i*dof+k];
1302:         }
1303:       } else {
1304:         for (k=0; k<dof; k++) {
1305:           mass_initial = mass_initial + hf*ptr_X0[i*dof+k];
1306:           mass_final = mass_final + hf*ptr_X[i*dof+k];
1307:         }
1308:       }
1309:     }
1310:     DMDAVecRestoreArrayRead(da,X0,(void*)&ptr_X0);
1311:     DMDAVecRestoreArrayRead(da,X,(void*)&ptr_X);
1312:     mass_difference = mass_final - mass_initial;
1313:     MPI_Allreduce(&mass_difference,&mass_differenceg,1,MPIU_SCALAR,MPIU_SUM,comm);
1314:     PetscPrintf(comm,"Mass difference %g\n",(double)mass_differenceg);
1315:     PetscPrintf(comm,"Final time %g, steps %D\n",(double)ptime,steps);
1316:     PetscPrintf(comm,"Maximum allowable stepsize according to CFL %g\n",(double)(1.0/ctx.cfl_idt));
1317:     if (ctx.exact) {
1318:       PetscReal nrm1=0;
1319:       SolutionErrorNorms_2WaySplit(&ctx,da,ptime,X,&nrm1);
1320:       PetscPrintf(comm,"Error ||x-x_e||_1 %g\n",(double)nrm1);
1321:     }
1322:     if (ctx.simulation) {
1323:       PetscReal    nrm1=0;
1324:       PetscViewer  fd;
1325:       char         filename[PETSC_MAX_PATH_LEN] = "binaryoutput";
1326:       Vec          XR;
1327:       PetscBool    flg;
1328:       const PetscScalar  *ptr_XR;
1329:       PetscOptionsGetString(NULL,NULL,"-f",filename,sizeof(filename),&flg);
1331:       PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_READ,&fd);
1332:       VecDuplicate(X0,&XR);
1333:       VecLoad(XR,fd);
1334:       PetscViewerDestroy(&fd);
1335:       VecGetArrayRead(X,&ptr_X);
1336:       VecGetArrayRead(XR,&ptr_XR);
1337:       for (i=xs;i<xs+xm;i++) {
1338:         if (i < ctx.sf || i > ctx.fs-1)
1339:           for (k=0; k<dof; k++) nrm1 = nrm1 + hs*PetscAbs(ptr_X[i*dof+k]-ptr_XR[i*dof+k]);
1340:         else
1341:           for (k=0; k<dof; k++) nrm1 = nrm1 + hf*PetscAbs(ptr_X[i*dof+k]-ptr_XR[i*dof+k]);
1342:       }
1343:       VecRestoreArrayRead(X,&ptr_X);
1344:       VecRestoreArrayRead(XR,&ptr_XR);
1345:       PetscPrintf(comm,"Error ||x-x_e||_1 %g\n",(double)nrm1);
1346:       VecDestroy(&XR);
1347:     }
1348:   }

1350:   SolutionStatsView(da,X,PETSC_VIEWER_STDOUT_WORLD);
1351:   if (draw & 0x1) VecView(X0,PETSC_VIEWER_DRAW_WORLD);
1352:   if (draw & 0x2) VecView(X,PETSC_VIEWER_DRAW_WORLD);
1353:   if (draw & 0x4) {
1354:     Vec Y;
1355:     VecDuplicate(X,&Y);
1356:     FVSample_2WaySplit(&ctx,da,ptime,Y);
1357:     VecAYPX(Y,-1,X);
1358:     VecView(Y,PETSC_VIEWER_DRAW_WORLD);
1359:     VecDestroy(&Y);
1360:   }

1362:   if (view_final) {
1363:     PetscViewer viewer;
1364:     PetscViewerASCIIOpen(PETSC_COMM_WORLD,final_fname,&viewer);
1365:     PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_MATLAB);
1366:     VecView(X,viewer);
1367:     PetscViewerPopFormat(viewer);
1368:     PetscViewerDestroy(&viewer);
1369:   }

1371:   /* Clean up */
1372:   (*ctx.physics2.destroy)(ctx.physics2.user);
1373:   for (i=0; i<ctx.physics2.dof; i++) PetscFree(ctx.physics2.fieldname[i]);
1374:   PetscFree(ctx.physics2.bcinflowindex);
1375:   PetscFree(ctx.ub);
1376:   PetscFree4(ctx.R,ctx.Rinv,ctx.cjmpLR,ctx.cslope);
1377:   PetscFree3(ctx.uLR,ctx.flux,ctx.speeds);
1378:   VecDestroy(&X);
1379:   VecDestroy(&X0);
1380:   VecDestroy(&R);
1381:   DMDestroy(&da);
1382:   TSDestroy(&ts);
1383:   ISDestroy(&ctx.iss);
1384:   ISDestroy(&ctx.isf);
1385:   ISDestroy(&ctx.issb);
1386:   PetscFree(index_slow);
1387:   PetscFree(index_fast);
1388:   PetscFree(index_slowbuffer);
1389:   PetscFunctionListDestroy(&limiters);
1390:   PetscFunctionListDestroy(&physics);
1391:   PetscFinalize();
1392:   return 0;
1393: }

1395: /*TEST

1397:     build:
1398:       requires: !complex !single
1399:       depends: finitevolume1d.c

1401:     test:
1402:       suffix: 1
1403:       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
1404:       output_file: output/ex4_1.out

1406:     test:
1407:       suffix: 2
1408:       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 1
1409:       output_file: output/ex4_1.out

1411:     test:
1412:       suffix: 3
1413:       args: -da_grid_x 40 -initial 1 -hratio 2 -limit mc -ts_dt 0.1 -ts_max_steps 24 -ts_max_time 7.0 -ts_type mprk -ts_mprk_type 2a22 -physics shallow -bc_type outflow -xmin 0 -xmax 50 -ts_use_splitrhsfunction 0
1414:       output_file: output/ex4_3.out

1416:     test:
1417:       suffix: 4
1418:       nsize: 2
1419:       args: -da_grid_x 40 -initial 1 -hratio 2 -limit mc -ts_dt 0.1 -ts_max_steps 24 -ts_max_time 7.0 -ts_type mprk -ts_mprk_type 2a22 -physics shallow -bc_type outflow -xmin 0 -xmax 50 -ts_use_splitrhsfunction 1
1420:       output_file: output/ex4_3.out

1422:     test:
1423:       suffix: 5
1424:       nsize: 4
1425:       args: args: -da_grid_x 40 -initial 1 -hratio 2 -limit mc -ts_dt 0.1 -ts_max_steps 24 -ts_max_time 7.0 -ts_type mprk -ts_mprk_type 2a22 -physics shallow -bc_type outflow -xmin 0 -xmax 50 -ts_use_splitrhsfunction 1
1426:       output_file: output/ex4_3.out
1427: TEST*/