Stuck Trying to Implement 3D Wave Equation in PyOpenCL

Question

I'm attempting to implement the discrete time wave equation in OpenCL. I think I'm pretty close, but the results look like what I would expect from the heat equation. I know they're very similar, but when I've implemented the 2D wave equation (not using OpenCL) I got distinct wavefronts and reflections. With the OpenCL kernel below everything diffuses until it is a wash.

 __kernel void wave_calc(
    __global float* height,
    __global float* height_old,
    const unsigned int len_x,
    const unsigned int len_y,
    const unsigned int len_z,
    const float dtxc_term)
{

    unsigned int x = get_global_id(0);
    unsigned int y = get_global_id(1);
    unsigned int z = get_global_id(2);
    int this_cell = x + len_y * (y + len_x * z);
    float laplacian;

    if (x==0 || x==(len_x-1) || y==0 || y==(len_y-1) || z==0 || z==(len_z-1)) {
       laplacian = 0;
       height_old[this_cell]  = height[this_cell];
       height[this_cell] = 0;
    }

    else if ( x < len_x-1 && y < len_y-1 && z < len_z-1 ){
        int n1 = x - 1 + len_y * (y + len_x * z);
        int n2 = x + 1 + len_y * (y + len_x * z);
        int n3 = x + len_y * (y - 1 + len_x * z);
        int n4 = x + len_y * (y + 1 + len_x * z);
        int n5 = x + len_y * (y + len_x * (z -1));
        int n6 = x + len_y * (y + len_x * (z + 1));

        laplacian = -6 * height[this_cell] +
                                    height[n1] +
                                    height[n2] +
                                    height[n3] +
                                    height[n4] +
                                    height[n5] +
                                    height[n6];

        height_old[this_cell]  = height[this_cell];
        height[this_cell] = (dtxc_term*laplacian+2*height[this_cell]) - height_old[this_cell];

    }  
}

(DTXC is the result of ((DT * DT)/(DX * DX)) * C passed from the host)

Every step I copy height back to the host for plotting, and then call the function again.

 for i in np.arange(steps):
        #copy height from host to device
        cl.enqueue_copy(queue, d_height, h_height)
        #step once
        wave_calc(queue, field_3d.shape, None, d_height, d_height_old, LEN_X, LEN_Y, LEN_Z, DTXC)

        queue.finish()
        #copy height back
        cl.enqueue_copy(queue, h_height, d_height)

#do my plotting

Any thoughts/suggestions/condescending remarks? All would be appreciated. :)

Here is an update to answer Joel's question:

I'm not much good when it comes to calculus, but I'm taking a working C++ implementation in 2D and trying to adapt it for 3D. Below is the C++. The only modification I made was to the loop, since there are 6 neighbor cells in 3D instead of 4. In both cases the outer walls of the plane/cube are set to 0:

for(int x=1; x<field.xRes()-1;x++) {
            for (int y=1; y<field.yRes()-1; y++) {

                    laplacian(x,y) = -4 * height(x,y) +
                    height(x-1,y) +
                    height(x+1,y) +
                    height(x,y-1) +
                    height(x,y+1);

            }
        }

const float dt = 0.001;
const float xLen = 1.0;
const float C = 1.0;
const float dx = xLen/xRes;

backup = height;
height = ((dt*dt)/(dx*dx))*C*laplacian+2*height;
height = height - heightOld;
heightOld = backup;