Actual source code: ex6.c

      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: */

      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";

     21: #include <petscts.h>
     22: #include <petscdm.h>
     23: #include <petscdmda.h>
     24: #include <petscdraw.h>
     25: #include "finitevolume1d.h"

     27: PETSC_STATIC_INLINE PetscReal RangeMod(PetscReal a,PetscReal xmin,PetscReal xmax) { PetscReal range = xmax-xmin; return xmin +PetscFmodReal(range+PetscFmodReal(a,range),range); }

     29: /* --------------------------------- Advection ----------------------------------- */
     30: typedef struct {
     31:   PetscReal a;                  /* advective velocity */
     32: } AdvectCtx;

     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;

     40:   speed     = ctx->a;
     41:   flux[0]   = PetscMax(0,speed)*uL[0] + PetscMin(0,speed)*uR[0];
     42:   *maxspeed = speed;
     43:   return(0);
     44: }

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

     51:   X[0]      = 1.;
     52:   Xi[0]     = 1.;
     53:   speeds[0] = ctx->a;
     54:   return(0);
     55: }

     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;

     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:   return(0);
     80: }

     82: static PetscErrorCode PhysicsCreate_Advect(FVCtx *ctx)
     83: {
     85:   AdvectCtx      *user;

     88:   PetscNew(&user);
     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:   PetscStrallocpy("u",&ctx->physics2.fieldname[0]);
     96:   user->a = 1;
     97:   PetscOptionsBegin(ctx->comm,ctx->prefix,"Options for advection","");
     98:   {
     99:     PetscOptionsReal("-physics_advect_a","Speed","",user->a,&user->a,NULL);
    100:   }
    101:   PetscOptionsEnd();
    102:   return(0);
    103: }

    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;

    114:   if (!ctx->physics2.sample2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Physics has not provided a sampling function");
    115:   DMDAGetInfo(da,0, &Mx,0,0, 0,0,0, &dof,0,0,0,0,0);
    116:   DMDAGetCorners(da,&xs,0,0,&xm,0,0);
    117:   DMDAVecGetArray(da,U,&u);
    118:   PetscMalloc1(dof,&uj);
    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:         (*ctx->physics2.sample2)(ctx->physics2.user,ctx->initial,ctx->bctype,ctx->xmin,ctx->xmax,time,xj,uj);
    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:         (*ctx->physics2.sample2)(ctx->physics2.user,ctx->initial,ctx->bctype,ctx->xmin,ctx->xmax,time,xj,uj);
    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:         (*ctx->physics2.sample2)(ctx->physics2.user,ctx->initial,ctx->bctype,ctx->xmin,ctx->xmax,time,xj,uj);
    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:   DMDAVecRestoreArray(da,U,&u);
    152:   PetscFree(uj);
    153:   return(0);
    154: }

    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;

    165:   VecGetSize(X,&Mx);
    166:   VecDuplicate(X,&Y);
    167:   FVSample_2WaySplit(ctx,da,t,Y);
    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:   VecGetArrayRead(X,&ptr_X);
    171:   VecGetArrayRead(Y,&ptr_Y);
    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:   VecRestoreArrayRead(X,&ptr_X);
    177:   VecRestoreArrayRead(Y,&ptr_Y);
    178:   VecDestroy(&Y);
    179:   return(0);
    180: }

    182: PetscErrorCode FVRHSFunction_2WaySplit(TS ts,PetscReal time,Vec X,Vec F,void *vctx)
    183: {
    184:   FVCtx          *ctx = (FVCtx*)vctx;
    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;

    193:   TSGetDM(ts,&da);
    194:   DMGetLocalVector(da,&Xloc);                          /* Xloc contains ghost points                                     */
    195:   DMDAGetInfo(da,0, &Mx,0,0, 0,0,0, &dof,0,0,0,0,0);   /* 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:   DMGlobalToLocalBegin(da,X,INSERT_VALUES,Xloc);       /* X is solution vector which does not contain ghost points       */
    199:   DMGlobalToLocalEnd(da,X,INSERT_VALUES,Xloc);

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

    203:   DMDAVecGetArray(da,Xloc,&x);
    204:   DMDAVecGetArray(da,F,&f);
    205:   DMDAGetArray(da,PETSC_TRUE,&slope);                  /* contains ghost points                                           */

    207:   DMDAGetCorners(da,&xs,0,0,&xm,0,0);

    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:     (*ctx->physics2.characteristic2)(ctx->physics2.user,dof,&x[i*dof],ctx->R,ctx->Rinv,ctx->speeds);
    222:     /* Evaluate jumps across interfaces (i-1, i) and (i, i+1), put in characteristic basis */
    223:     PetscArrayzero(ctx->cjmpLR,2*dof);
    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:   }

    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:       (*ctx->physics2.riemann2)(ctx->physics2.user,dof,uL,uR,ctx->flux,&maxspeed);
    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:       (*ctx->physics2.riemann2)(ctx->physics2.user,dof,uL,uR,ctx->flux,&maxspeed);
    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:       (*ctx->physics2.riemann2)(ctx->physics2.user,dof,uL,uR,ctx->flux,&maxspeed);
    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:       (*ctx->physics2.riemann2)(ctx->physics2.user,dof,uL,uR,ctx->flux,&maxspeed);
    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:       (*ctx->physics2.riemann2)(ctx->physics2.user,dof,uL,uR,ctx->flux,&maxspeed);
    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:   DMDAVecRestoreArray(da,Xloc,&x);
    316:   DMDAVecRestoreArray(da,F,&f);
    317:   DMDARestoreArray(da,PETSC_TRUE,&slope);
    318:   DMRestoreLocalVector(da,&Xloc);
    319:   MPI_Allreduce(&cfl_idt,&ctx->cfl_idt,1,MPIU_SCALAR,MPIU_MAX,PetscObjectComm((PetscObject)da));
    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:     TSGetTimeStep(ts,&dt);
    324:     TSGetTime(ts,&tnow);
    325:     if (dt > 0.5/ctx->cfl_idt) {
    326:       if (1) {
    327:         PetscPrintf(ctx->comm,"Stability constraint exceeded at t=%g, dt %g > %g\n",(double)tnow,(double)dt,(double)(0.5/ctx->cfl_idt));
    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:   return(0);
    332: }

    334: /* --------------------------------- Finite Volume Solver for slow components ----------------------------------- */
    335: PetscErrorCode FVRHSFunctionslow_2WaySplit(TS ts,PetscReal time,Vec X,Vec F,void *vctx)
    336: {
    337:   FVCtx          *ctx = (FVCtx*)vctx;
    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;

    346:   TSGetDM(ts,&da);
    347:   DMGetLocalVector(da,&Xloc);
    348:   DMDAGetInfo(da,0, &Mx,0,0, 0,0,0, &dof,0,0,0,0,0);
    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:   DMGlobalToLocalBegin(da,X,INSERT_VALUES,Xloc);
    352:   DMGlobalToLocalEnd(da,X,INSERT_VALUES,Xloc);
    353:   VecZeroEntries(F);
    354:   DMDAVecGetArray(da,Xloc,&x);
    355:   VecGetArray(F,&f);
    356:   DMDAGetArray(da,PETSC_TRUE,&slope);
    357:   DMDAGetCorners(da,&xs,0,0,&xm,0,0);

    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:       (*ctx->physics2.characteristic2)(ctx->physics2.user,dof,&x[i*dof],ctx->R,ctx->Rinv,ctx->speeds);
    373:       /* Evaluate jumps across interfaces (i-1, i) and (i, i+1), put in characteristic basis */
    374:       PetscArrayzero(ctx->cjmpLR,2*dof);
    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:   }

    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:       (*ctx->physics2.riemann2)(ctx->physics2.user,dof,uL,uR,ctx->flux,&maxspeed);
    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:       (*ctx->physics2.riemann2)(ctx->physics2.user,dof,uL,uR,ctx->flux,&maxspeed);
    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:       (*ctx->physics2.riemann2)(ctx->physics2.user,dof,uL,uR,ctx->flux,&maxspeed);
    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:       (*ctx->physics2.riemann2)(ctx->physics2.user,dof,uL,uR,ctx->flux,&maxspeed);
    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:   DMDAVecRestoreArray(da,Xloc,&x);
    456:   VecRestoreArray(F,&f);
    457:   DMDARestoreArray(da,PETSC_TRUE,&slope);
    458:   DMRestoreLocalVector(da,&Xloc);
    459:   MPI_Allreduce(&cfl_idt,&ctx->cfl_idt,1,MPIU_SCALAR,MPIU_MAX,PetscObjectComm((PetscObject)da));
    460:   return(0);
    461: }

    463: PetscErrorCode FVRHSFunctionslowbuffer_2WaySplit(TS ts,PetscReal time,Vec X,Vec F,void *vctx)
    464: {
    465:   FVCtx          *ctx = (FVCtx*)vctx;
    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;

    474:   TSGetDM(ts,&da);
    475:   DMGetLocalVector(da,&Xloc);
    476:   DMDAGetInfo(da,0, &Mx,0,0, 0,0,0, &dof,0,0,0,0,0);
    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:   DMGlobalToLocalBegin(da,X,INSERT_VALUES,Xloc);
    480:   DMGlobalToLocalEnd(da,X,INSERT_VALUES,Xloc);
    481:   VecZeroEntries(F);
    482:   DMDAVecGetArray(da,Xloc,&x);
    483:   VecGetArray(F,&f);
    484:   DMDAGetArray(da,PETSC_TRUE,&slope);
    485:   DMDAGetCorners(da,&xs,0,0,&xm,0,0);

    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:       (*ctx->physics2.characteristic2)(ctx->physics2.user,dof,&x[i*dof],ctx->R,ctx->Rinv,ctx->speeds);
    501:       /* Evaluate jumps across interfaces (i-1, i) and (i, i+1), put in characteristic basis */
    502:       PetscArrayzero(ctx->cjmpLR,2*dof);
    503:       cjmpL = &ctx->cjmpLR[0];
    504:       cjmpR = &ctx->cjmpLR[dof];
    505:       for (j=0; j<dof; j++) {
    506:         PetscScalar jmpL,jmpR;
    507:         jmpL = x[(i+0)*dof+j]-x[(i-1)*dof+j];
    508:         jmpR = x[(i+1)*dof+j]-x[(i+0)*dof+j];
    509:         for (k=0; k<dof; k++) {
    510:           cjmpL[k] += ctx->Rinv[k+j*dof]*jmpL;
    511:           cjmpR[k] += ctx->Rinv[k+j*dof]*jmpR;
    512:         }
    513:       }
    514:       /* Apply limiter to the left and right characteristic jumps */
    515:       info.m  = dof;
    516:       info.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:   }

    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:       (*ctx->physics2.riemann2)(ctx->physics2.user,dof,uL,uR,ctx->flux,&maxspeed);
    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:       (*ctx->physics2.riemann2)(ctx->physics2.user,dof,uL,uR,ctx->flux,&maxspeed);
    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:       (*ctx->physics2.riemann2)(ctx->physics2.user,dof,uL,uR,ctx->flux,&maxspeed);
    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:       (*ctx->physics2.riemann2)(ctx->physics2.user,dof,uL,uR,ctx->flux,&maxspeed);
    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:       (*ctx->physics2.riemann2)(ctx->physics2.user,dof,uL,uR,ctx->flux,&maxspeed);
    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:       (*ctx->physics2.riemann2)(ctx->physics2.user,dof,uL,uR,ctx->flux,&maxspeed);
    598:       if (i > xs) {
    599:         for (j=0; j<dof; j++) f[(islow-1)*dof+j] -= ctx->flux[j]/hxs;
    600:       }
    601:     }
    602:   }
    603:   DMDAVecRestoreArray(da,Xloc,&x);
    604:   VecRestoreArray(F,&f);
    605:   DMDARestoreArray(da,PETSC_TRUE,&slope);
    606:   DMRestoreLocalVector(da,&Xloc);
    607:   return(0);
    608: }

    610: /* --------------------------------- Finite Volume Solver for fast  parts ----------------------------------- */
    611: PetscErrorCode FVRHSFunctionfast_2WaySplit(TS ts,PetscReal time,Vec X,Vec F,void *vctx)
    612: {
    613:   FVCtx          *ctx = (FVCtx*)vctx;
    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;

    622:   TSGetDM(ts,&da);
    623:   DMGetLocalVector(da,&Xloc);
    624:   DMDAGetInfo(da,0, &Mx,0,0, 0,0,0, &dof,0,0,0,0,0);
    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:   DMGlobalToLocalBegin(da,X,INSERT_VALUES,Xloc);
    628:   DMGlobalToLocalEnd(da,X,INSERT_VALUES,Xloc);
    629:   VecZeroEntries(F);
    630:   DMDAVecGetArray(da,Xloc,&x);
    631:   VecGetArray(F,&f);
    632:   DMDAGetArray(da,PETSC_TRUE,&slope);
    633:   DMDAGetCorners(da,&xs,0,0,&xm,0,0);

    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:       (*ctx->physics2.characteristic2)(ctx->physics2.user,dof,&x[i*dof],ctx->R,ctx->Rinv,ctx->speeds);
    648:       PetscArrayzero(ctx->cjmpLR,2*dof);
    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:   }

    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:       (*ctx->physics2.riemann2)(ctx->physics2.user,dof,uL,uR,ctx->flux,&maxspeed);
    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:       (*ctx->physics2.riemann2)(ctx->physics2.user,dof,uL,uR,ctx->flux,&maxspeed);
    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:       (*ctx->physics2.riemann2)(ctx->physics2.user,dof,uL,uR,ctx->flux,&maxspeed);
    709:       if (i > xs) {
    710:         for (j=0; j<dof; j++) f[(ifast-1)*dof+j] -= ctx->flux[j]/hxf;
    711:       }
    712:     }
    713:   }
    714:   DMDAVecRestoreArray(da,Xloc,&x);
    715:   VecRestoreArray(F,&f);
    716:   DMDARestoreArray(da,PETSC_TRUE,&slope);
    717:   DMRestoreLocalVector(da,&Xloc);
    718:   return(0);
    719: }

    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;

    735:   PetscInitialize(&argc,&argv,0,help);if (ierr) return ierr;
    736:   comm = PETSC_COMM_WORLD;
    737:   PetscMemzero(&ctx,sizeof(ctx));

    739:   /* Register limiters to be available on the command line */
    740:   PetscFunctionListAdd(&limiters,"upwind"              ,Limit2_Upwind);
    741:   PetscFunctionListAdd(&limiters,"lax-wendroff"        ,Limit2_LaxWendroff);
    742:   PetscFunctionListAdd(&limiters,"beam-warming"        ,Limit2_BeamWarming);
    743:   PetscFunctionListAdd(&limiters,"fromm"               ,Limit2_Fromm);
    744:   PetscFunctionListAdd(&limiters,"minmod"              ,Limit2_Minmod);
    745:   PetscFunctionListAdd(&limiters,"superbee"            ,Limit2_Superbee);
    746:   PetscFunctionListAdd(&limiters,"mc"                  ,Limit2_MC);
    747:   PetscFunctionListAdd(&limiters,"koren3"              ,Limit2_Koren3);

    749:   /* Register physical models to be available on the command line */
    750:   PetscFunctionListAdd(&physics,"advect"          ,PhysicsCreate_Advect);

    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:   PetscOptionsBegin(comm,NULL,"Finite Volume solver options","");
    758:   PetscOptionsReal("-xmin","X min","",ctx.xmin,&ctx.xmin,NULL);
    759:   PetscOptionsReal("-xmax","X max","",ctx.xmax,&ctx.xmax,NULL);
    760:   PetscOptionsFList("-limit","Name of flux imiter to use","",limiters,lname,lname,sizeof(lname),NULL);
    761:   PetscOptionsInt("-draw","Draw solution vector, bitwise OR of (1=initial,2=final,4=final error)","",draw,&draw,NULL);
    762:   PetscOptionsString("-view_final","Write final solution in ASCII MATLAB format to given file name","",final_fname,final_fname,sizeof(final_fname),&view_final);
    763:   PetscOptionsInt("-initial","Initial condition (depends on the physics)","",ctx.initial,&ctx.initial,NULL);
    764:   PetscOptionsBool("-exact","Compare errors with exact solution","",ctx.exact,&ctx.exact,NULL);
    765:   PetscOptionsBool("-simulation","Compare errors with reference solution","",ctx.simulation,&ctx.simulation,NULL);
    766:   PetscOptionsReal("-cfl","CFL number to time step at","",ctx.cfl,&ctx.cfl,NULL);
    767:   PetscOptionsEnum("-bc_type","Boundary condition","",FVBCTypes,(PetscEnum)ctx.bctype,(PetscEnum*)&ctx.bctype,NULL);
    768:   PetscOptionsInt("-hratio","Spacing ratio","",ctx.hratio,&ctx.hratio,NULL);
    769:   PetscOptionsEnd();

    771:   /* Choose the limiter from the list of registered limiters */
    772:   PetscFunctionListFind(limiters,lname,&ctx.limit2);
    773:   if (!ctx.limit2) SETERRQ1(PETSC_COMM_SELF,1,"Limiter '%s' not found",lname);

    775:   /* Choose the physics from the list of registered models */
    776:   {
    777:     PetscErrorCode (*r)(FVCtx*);
    778:     PetscFunctionListFind(physics,physname,&r);
    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:     (*r)(&ctx);
    782:   }

    784:   /* Create a DMDA to manage the parallel grid */
    785:   DMDACreate1d(comm,DM_BOUNDARY_PERIODIC,50,ctx.physics2.dof,2,NULL,&da);
    786:   DMSetFromOptions(da);
    787:   DMSetUp(da);
    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:     DMDASetFieldName(da,i,ctx.physics2.fieldname[i]);
    792:   }
    793:   DMDAGetInfo(da,0, &Mx,0,0, 0,0,0, &dof,0,0,0,0,0);
    794:   DMDAGetCorners(da,&xs,0,0,&xm,0,0);

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

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

    803:   /* Create a vector to store the solution and to save the initial state */
    804:   DMCreateGlobalVector(da,&X);
    805:   VecDuplicate(X,&X0);
    806:   VecDuplicate(X,&R);

    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:   PetscMalloc1(xm*dof,&index_slow);
    816:   PetscMalloc1(xm*dof,&index_fast);
    817:   PetscMalloc1(6*dof,&index_slowbuffer);
    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:   ISCreateGeneral(PETSC_COMM_WORLD,islow,index_slow,PETSC_COPY_VALUES,&ctx.iss);
    834:   ISCreateGeneral(PETSC_COMM_WORLD,ifast,index_fast,PETSC_COPY_VALUES,&ctx.isf);
    835:   ISCreateGeneral(PETSC_COMM_WORLD,islowbuffer,index_slowbuffer,PETSC_COPY_VALUES,&ctx.issb);

    837:   /* Create a time-stepping object */
    838:   TSCreate(comm,&ts);
    839:   TSSetDM(ts,da);
    840:   TSSetRHSFunction(ts,R,FVRHSFunction_2WaySplit,&ctx);
    841:   TSRHSSplitSetIS(ts,"slow",ctx.iss);
    842:   TSRHSSplitSetIS(ts,"slowbuffer",ctx.issb);
    843:   TSRHSSplitSetIS(ts,"fast",ctx.isf);
    844:   TSRHSSplitSetRHSFunction(ts,"slow",NULL,FVRHSFunctionslow_2WaySplit,&ctx);
    845:   TSRHSSplitSetRHSFunction(ts,"fast",NULL,FVRHSFunctionfast_2WaySplit,&ctx);
    846:   TSRHSSplitSetRHSFunction(ts,"slowbuffer",NULL,FVRHSFunctionslowbuffer_2WaySplit,&ctx);

    848:   TSSetType(ts,TSSSP);
    849:   /*TSSetType(ts,TSMPRK);*/
    850:   TSSetMaxTime(ts,10);
    851:   TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);

    853:   /* Compute initial conditions and starting time step */
    854:   FVSample_2WaySplit(&ctx,da,0,X0);
    855:   FVRHSFunction_2WaySplit(ts,0,X0,X,(void*)&ctx); /* Initial function evaluation, only used to determine max speed */
    856:   VecCopy(X0,X);                        /* The function value was not used so we set X=X0 again */
    857:   TSSetTimeStep(ts,ctx.cfl/ctx.cfl_idt);
    858:   TSSetFromOptions(ts); /* Take runtime options */
    859:   SolutionStatsView(da,X,PETSC_VIEWER_STDOUT_WORLD);
    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;

    867:     TSSolve(ts,X);
    868:     TSGetSolveTime(ts,&ptime);
    869:     TSGetStepNumber(ts,&steps);
    870:     /* calculate the total mass at initial time and final time */
    871:     mass_initial = 0.0;
    872:     mass_final   = 0.0;
    873:     DMDAVecGetArrayRead(da,X0,(void*)&ptr_X0);
    874:     DMDAVecGetArrayRead(da,X,(void*)&ptr_X);
    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:     DMDAVecRestoreArrayRead(da,X0,(void*)&ptr_X0);
    889:     DMDAVecRestoreArrayRead(da,X,(void*)&ptr_X);
    890:     mass_difference = mass_final - mass_initial;
    891:     MPI_Allreduce(&mass_difference,&mass_differenceg,1,MPIU_SCALAR,MPIU_SUM,comm);
    892:     PetscPrintf(comm,"Mass difference %g\n",(double)mass_differenceg);
    893:     PetscPrintf(comm,"Final time %g, steps %D\n",(double)ptime,steps);
    894:     PetscPrintf(comm,"Maximum allowable stepsize according to CFL %g\n",(double)(1.0/ctx.cfl_idt));
    895:     if (ctx.exact) {
    896:       PetscReal nrm1=0;
    897:       SolutionErrorNorms_2WaySplit(&ctx,da,ptime,X,&nrm1);
    898:       PetscPrintf(comm,"Error ||x-x_e||_1 %g\n",(double)nrm1);
    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:       PetscOptionsGetString(NULL,NULL,"-f",filename,sizeof(filename),&flg);
    908:       if (!flg) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_USER,"Must indicate binary file with the -f option");
    909:       PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_READ,&fd);
    910:       VecDuplicate(X0,&XR);
    911:       VecLoad(XR,fd);
    912:       PetscViewerDestroy(&fd);
    913:       VecGetArrayRead(X,&ptr_X);
    914:       VecGetArrayRead(XR,&ptr_XR);
    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:       VecRestoreArrayRead(X,&ptr_X);
    922:       VecRestoreArrayRead(XR,&ptr_XR);
    923:       PetscPrintf(comm,"Error ||x-x_e||_1 %g\n",(double)nrm1);
    924:       VecDestroy(&XR);
    925:     }
    926:   }

    928:   SolutionStatsView(da,X,PETSC_VIEWER_STDOUT_WORLD);
    929:   if (draw & 0x1) {VecView(X0,PETSC_VIEWER_DRAW_WORLD);}
    930:   if (draw & 0x2) {VecView(X,PETSC_VIEWER_DRAW_WORLD);}
    931:   if (draw & 0x4) {
    932:     Vec Y;
    933:     VecDuplicate(X,&Y);
    934:     FVSample_2WaySplit(&ctx,da,ptime,Y);
    935:     VecAYPX(Y,-1,X);
    936:     VecView(Y,PETSC_VIEWER_DRAW_WORLD);
    937:     VecDestroy(&Y);
    938:   }

    940:   if (view_final) {
    941:     PetscViewer viewer;
    942:     PetscViewerASCIIOpen(PETSC_COMM_WORLD,final_fname,&viewer);
    943:     PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_MATLAB);
    944:     VecView(X,viewer);
    945:     PetscViewerPopFormat(viewer);
    946:     PetscViewerDestroy(&viewer);
    947:   }

    949:   /* Clean up */
    950:   (*ctx.physics2.destroy)(ctx.physics2.user);
    951:   for (i=0; i<ctx.physics2.dof; i++) {PetscFree(ctx.physics2.fieldname[i]);}
    952:   PetscFree4(ctx.R,ctx.Rinv,ctx.cjmpLR,ctx.cslope);
    953:   PetscFree3(ctx.uLR,ctx.flux,ctx.speeds);
    954:   VecDestroy(&X);
    955:   VecDestroy(&X0);
    956:   VecDestroy(&R);
    957:   DMDestroy(&da);
    958:   TSDestroy(&ts);
    959:   ISDestroy(&ctx.iss);
    960:   ISDestroy(&ctx.isf);
    961:   ISDestroy(&ctx.issb);
    962:   PetscFree(index_slow);
    963:   PetscFree(index_fast);
    964:   PetscFree(index_slowbuffer);
    965:   PetscFunctionListDestroy(&limiters);
    966:   PetscFunctionListDestroy(&physics);
    967:   PetscFinalize();
    968:   return ierr;
    969: }

    971: /*TEST

    973:     build:
    974:       requires: !complex
    975:       depends: finitevolume1d.c

    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

    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

    986: TEST*/