     // Get X componet of gravity field
     id_l = (xid-1) + (yid)*nx + (zid)*nx*ny;
     id_r = (xid+1) + (yid)*nx + (zid)*nx*ny;
-    // pot_l = dev_conserved[field_pot*n_cells + id_l];
-    // pot_r = dev_conserved[field_pot*n_cells + id_r];
     pot_l = dev_potential[id_l];
     pot_r = dev_potential[id_r];
     gx = -0.5*( pot_r - pot_l ) / dx;

    //Get Y componet of gravity field
    id_l = (xid) + (yid-1)*nx + (zid)*nx*ny;
    id_r = (xid) + (yid+1)*nx + (zid)*nx*ny;
    // pot_l = dev_conserved[field_pot*n_cells + id_l];
    // pot_r = dev_conserved[field_pot*n_cells + id_r];
    pot_l = dev_potential[id_l];
    pot_r = dev_potential[id_r];
    gy = -0.5*( pot_r - pot_l ) / dy;

    //Get Z componet of gravity field
    id_l = (xid) + (yid)*nx + (zid-1)*nx*ny;
    id_r = (xid) + (yid)*nx + (zid+1)*nx*ny;
    pot_l = dev_potential[id_l];
    pot_r = dev_potential[id_r];
    gz = -0.5*( pot_r - pot_l ) / dz;
