     vx_i =  dev_conserved[o1*n_cells + id] / d_i;
     vy_i =  dev_conserved[o2*n_cells + id] / d_i;
     vz_i =  dev_conserved[o3*n_cells + id] / d_i;
+    #ifdef DE //PRESSURE_DE
+    E = dev_conserved[4*n_cells + id];
+    GE = dev_conserved[(n_fields-1)*n_cells + id];
+    E_kin = 0.5 * d_i * ( vx_i*vx_i + vy_i*vy_i + vz_i*vz_i );
+    p_i = Get_Pressure_From_DE( E, E - E_kin, GE, gamma ); 
+    #else
     p_i  = (dev_conserved[4*n_cells + id] - 0.5*d_i*(vx_i*vx_i + vy_i*vy_i + vz_i*vz_i)) * (gamma - 1.0);
+    #endif
     p_i  = fmax(p_i, (Real) TINY_NUMBER);
     #ifdef DE
     ge_i = dev_conserved[(n_fields-1)*n_cells + id] / d_i;

    vx_ipt =  dev_conserved[o1*n_cells + id] / d_ipt;
    vy_ipt =  dev_conserved[o2*n_cells + id] / d_ipt;
    vz_ipt =  dev_conserved[o3*n_cells + id] / d_ipt;
    p_ipt  = (dev_conserved[4*n_cells + id] - 0.5*d_ipt*(vx_ipt*vx_ipt + vy_ipt*vy_ipt + vz_ipt*vz_ipt)) * (gamma - 1.0);
    p_ipt  = fmax(p_ipt, (Real) TINY_NUMBER);
    #ifdef DE
    ge_ipt = dev_conserved[(n_fields-1)*n_cells + id] / d_ipt;

  int xid = tid - zid*nx*ny - yid*nx;

  int xs, xe, ys, ye, zs, ze;
  // if (dir == 0) {
  //   xs = 3; xe = nx-4;
  //   ys = 0; ye = ny;
  //   zs = 0; ze = nz;
  // }
  // if (dir == 1) {
  //   xs = 0; xe = nx;
  //   ys = 3; ye = ny-4;
  //   zs = 0; ze = nz;
  // }
  // if (dir == 2) {
  //   xs = 0; xe = nx;
  //   ys = 0; ye = ny;
  //   zs = 3; ze = nz-4;
  // }
  // NOTE: There are 2 ghost cells, not 3. Flattening not allowed!
  if (dir == 0) {
    xs = 2; xe = nx-3;
    ys = 0; ye = ny;
    zs = 0; ze = nz;
  }
  if (dir == 1) {
    xs = 0; xe = nx;
    ys = 2; ye = ny-3;
    zs = 0; ze = nz;
  }
  if (dir == 2) {
    xs = 0; xe = nx;
    ys = 0; ye = ny;
    zs = 2; ze = nz-3;
  }

  if (xid >= xs && xid < xe && yid >= ys && yid < ye && zid >= zs && zid < ze)
