creating branch for adding support for CUDA 3.x and Fermi

[git-p4: depot-paths = "//depot/gpgpu_sim_research/fermi/distribution/": change = 6829]

1 files changed, 551 insertions, 0 deletions

diff --git a/benchmarks/CUDA/DG/src/MaxwellsKernel3d.cu b/benchmarks/CUDA/DG/src/MaxwellsKernel3d.cu
new file mode 100644
index 0000000..bfa4396
--- /dev/null
+++ b/benchmarks/CUDA/DG/src/MaxwellsKernel3d.cu
@@ -0,0 +1,551 @@
+/* -*- mode: C; c-basic-offset: 8; c-indent-level: 8; c-continued-statement-offset: 8; c-label-offset: -8; -*- */
+
+#include <stdio.h>
+#include <cuda.h>
+
+texture<float4, 1, cudaReadModeElementType> t_LIFT;
+texture<float4, 1, cudaReadModeElementType> t_DrDsDt;
+texture<float, 1, cudaReadModeElementType> t_Dr;
+texture<float, 1, cudaReadModeElementType> t_Ds;
+texture<float, 1, cudaReadModeElementType> t_Dt;
+texture<float, 1, cudaReadModeElementType> t_vgeo;
+texture<float4, 1, cudaReadModeElementType> t_vgeo4;
+texture<float, 1, cudaReadModeElementType> t_Q;
+texture<float, 1, cudaReadModeElementType> t_partQ;
+texture<float, 1, cudaReadModeElementType> t_surfinfo;
+
+static float *c_LIFT;
+static float *c_DrDsDt;
+static float *c_surfinfo;
+static float *c_vgeo;
+static float *c_Q; 
+static float *c_partQ; 
+static float *c_rhsQ; 
+static float *c_resQ; 
+static float *c_tmp;
+
+extern "C"
+{
+
+#include "fem.h"
+
+double InitGPU3d(Mesh *mesh, int Nfields){
+
+  /* Q  */
+  int sz = mesh->K*(BSIZE)*p_Nfields*sizeof(float); 
+
+  float *f_Q = (float*) calloc(mesh->K*BSIZE*p_Nfields, sizeof(float));
+  cudaMalloc  ((void**) &c_Q, sz);
+  cudaMalloc  ((void**) &c_rhsQ, sz);
+  cudaMalloc  ((void**) &c_resQ, sz);
+  cudaMalloc  ((void**) &c_tmp, sz);
+  cudaMemcpy( c_Q,    f_Q, sz, cudaMemcpyHostToDevice);
+  cudaMemcpy( c_rhsQ, f_Q, sz, cudaMemcpyHostToDevice);
+  cudaMemcpy( c_resQ, f_Q, sz, cudaMemcpyHostToDevice);
+  cudaMemcpy( c_tmp,  f_Q, sz, cudaMemcpyHostToDevice);
+
+  cudaBindTexture(0,  t_Q, c_Q, sz);
+
+  sz = mesh->parNtotalout*sizeof(float);
+  cudaMalloc((void**) &c_partQ, sz);
+  cudaBindTexture(0,  t_partQ, c_partQ, sz);
+
+  /*  LIFT  */
+   sz = p_Np*(p_Nfp)*p_Nfaces*sizeof(float);
+#if 0
+   float *f_LIFT = (float*) malloc(sz);
+   int skL = 0;
+   for(int m=0;m<p_Nfp*p_Nfaces;++m){
+     for(int n=0;n<p_Np;++n){
+       f_LIFT[skL++] = d_LIFT[n+p_Np*m];
+     }
+   }
+#else
+   float *f_LIFT = (float*) malloc(sz);
+   int skL = 0;
+   for(int m=0;m<p_Nfp;++m){
+     for(int n=0;n<p_Np;++n){
+       for(int f=0;f<p_Nfaces;++f){
+	 f_LIFT[skL++] = mesh->LIFT[0][p_Nfp*p_Nfaces*n+(f+p_Nfaces*m)];
+       }
+     }
+   }
+#endif
+   cudaMalloc  ((void**) &c_LIFT, sz);
+   cudaMemcpy( c_LIFT, f_LIFT, sz, cudaMemcpyHostToDevice);
+   
+   /* Bind the array to the texture */
+   cudaBindTexture(0,  t_LIFT, c_LIFT, sz);
+
+   /* DrDsDt */
+   sz = BSIZE*BSIZE*4*sizeof(float);
+
+   float* h_DrDsDt = (float*) calloc(BSIZE*BSIZE, sizeof(float4));
+   int sk = 0;
+   /* note transposed arrays to avoid "bank conflicts" */
+   for(int n=0;n<p_Np;++n){
+     for(int m=0;m<p_Np;++m){
+       h_DrDsDt[4*(m+n*BSIZE)+0] = mesh->Dr[0][n+m*p_Np];
+       h_DrDsDt[4*(m+n*BSIZE)+1] = mesh->Ds[0][n+m*p_Np];
+       h_DrDsDt[4*(m+n*BSIZE)+2] = mesh->Dt[0][n+m*p_Np];
+     }
+   }
+	   
+   cudaMalloc  ((void**) &c_DrDsDt, sz);
+   cudaMemcpy( c_DrDsDt, h_DrDsDt, sz, cudaMemcpyHostToDevice);
+   
+   /* Bind the array to the texture */
+   cudaBindTexture(0,  t_DrDsDt, c_DrDsDt, sz);
+
+   free(h_DrDsDt);
+
+   /* vgeo */
+   double drdx, dsdx, dtdx;
+   double drdy, dsdy, dtdy;
+   double drdz, dsdz, dtdz, J;
+   float *vgeo = (float*) calloc(12*mesh->K, sizeof(float));
+
+   for(int k=0;k<mesh->K;++k){
+     GeometricFactors3d(mesh, k, 
+			&drdx, &dsdx, &dtdx,
+			&drdy, &dsdy, &dtdy,
+			&drdz, &dsdz, &dtdz, &J);
+
+     vgeo[k*12+0] = drdx; vgeo[k*12+1] = drdy; vgeo[k*12+2] = drdz;
+     vgeo[k*12+4] = dsdx; vgeo[k*12+5] = dsdy; vgeo[k*12+6] = dsdz;
+     vgeo[k*12+8] = dtdx; vgeo[k*12+9] = dtdy; vgeo[k*12+10] = dtdz;
+
+   }
+
+   sz = mesh->K*12*sizeof(float);
+   cudaMalloc  ((void**) &c_vgeo, sz);
+   cudaMemcpy( c_vgeo, vgeo, sz, cudaMemcpyHostToDevice);
+   cudaBindTexture(0,  t_vgeo, c_vgeo, sz);
+   
+   /* surfinfo (vmapM, vmapP, Fscale, Bscale, nx, ny, nz, 0) */
+   sz = mesh->K*p_Nfp*p_Nfaces*7*sizeof(float); 
+   float* h_surfinfo = (float*) malloc(sz); 
+   
+   /* local-local info */
+   sk = 0;
+   int skP = -1;
+   double *nxk = BuildVector(mesh->Nfaces);
+   double *nyk = BuildVector(mesh->Nfaces);
+   double *nzk = BuildVector(mesh->Nfaces);
+   double *sJk = BuildVector(mesh->Nfaces);
+
+   double dt = 1e6;
+
+   for(int k=0;k<mesh->K;++k){
+
+     GeometricFactors3d(mesh, k, 
+			&drdx, &dsdx, &dtdx,
+			&drdy, &dsdy, &dtdy,
+			&drdz, &dsdz, &dtdz, &J);
+
+     Normals3d(mesh, k, nxk, nyk, nzk, sJk);
+     
+     for(int f=0;f<mesh->Nfaces;++f){
+
+	     dt = min(dt, J/sJk[f]);
+       
+       for(int m=0;m<p_Nfp;++m){
+	 int n = m + f*p_Nfp + p_Nfp*p_Nfaces*k;
+	 int idM = mesh->vmapM[n];
+	 int idP = mesh->vmapP[n];
+	 int  nM = idM%p_Np; 
+	 int  nP = idP%p_Np; 
+	 int  kM = (idM-nM)/p_Np;
+	 int  kP = (idP-nP)/p_Np;
+	 idM = nM + Nfields*BSIZE*kM;
+	 idP = nP + Nfields*BSIZE*kP;
+	 
+	 /* stub resolve some other way */
+	 if(mesh->vmapP[n]<0){
+	   idP = mesh->vmapP[n]; /* -ve numbers */
+	 }
+ 
+	 sk = 7*p_Nfp*p_Nfaces*k+m+f*p_Nfp;
+	 h_surfinfo[sk + 0*p_Nfp*p_Nfaces] = idM;
+	 h_surfinfo[sk + 1*p_Nfp*p_Nfaces] = idP;
+	 h_surfinfo[sk + 2*p_Nfp*p_Nfaces] = sJk[f]/(2.*J);
+	 h_surfinfo[sk + 3*p_Nfp*p_Nfaces] = (idM==idP)?-1.:1.;
+	 h_surfinfo[sk + 4*p_Nfp*p_Nfaces] = nxk[f];
+	 h_surfinfo[sk + 5*p_Nfp*p_Nfaces] = nyk[f];
+	 h_surfinfo[sk + 6*p_Nfp*p_Nfaces] = nzk[f];
+       }
+     }
+   }
+	   
+   cudaMalloc  ((void**) &c_surfinfo, sz);
+   cudaMemcpy( c_surfinfo, h_surfinfo, sz, cudaMemcpyHostToDevice);
+
+   cudaBindTexture(0,  t_surfinfo, c_surfinfo, sz);
+
+   free(h_surfinfo);
+
+   sz = mesh->parNtotalout*sizeof(int);
+   cudaMalloc((void**) &(mesh->c_parmapOUT), sz);
+   cudaMemcpy(mesh->c_parmapOUT,  mesh->parmapOUT, sz, cudaMemcpyHostToDevice);
+
+   return dt;
+}
+
+
+
+__global__ void MaxwellsGPU_VOL_Kernel3D(float *g_rhsQ){
+
+  /* fastest */
+  __device__ __shared__ float s_Q[p_Nfields*BSIZE];
+  __device__ __shared__ float s_facs[12];
+
+  const int n = threadIdx.x;
+  const int k = blockIdx.x;
+  
+  /* "coalesced"  */
+  int m = n+k*p_Nfields*BSIZE;
+  int id = n;
+  s_Q[id] = tex1Dfetch(t_Q, m); m+=BSIZE; id+=BSIZE;
+  s_Q[id] = tex1Dfetch(t_Q, m); m+=BSIZE; id+=BSIZE;
+  s_Q[id] = tex1Dfetch(t_Q, m); m+=BSIZE; id+=BSIZE;
+  s_Q[id] = tex1Dfetch(t_Q, m); m+=BSIZE; id+=BSIZE;
+  s_Q[id] = tex1Dfetch(t_Q, m); m+=BSIZE; id+=BSIZE;
+  s_Q[id] = tex1Dfetch(t_Q, m); 
+
+  if(p_Np<12 && n==0)
+    for(m=0;m<12;++m)
+      s_facs[m] = tex1Dfetch(t_vgeo, 12*k+m);
+  else if(n<12 && p_Np>=12)
+    s_facs[n] = tex1Dfetch(t_vgeo, 12*k+n);
+
+  __syncthreads();
+
+  float dHxdr=0,dHxds=0,dHxdt=0;
+  float dHydr=0,dHyds=0,dHydt=0;
+  float dHzdr=0,dHzds=0,dHzdt=0;
+  float dExdr=0,dExds=0,dExdt=0;
+  float dEydr=0,dEyds=0,dEydt=0;
+  float dEzdr=0,dEzds=0,dEzdt=0;
+  float Q;
+
+  for(m=0;p_Np-m;){
+    float4 D = tex1Dfetch(t_DrDsDt, n+m*BSIZE);
+
+    id = m;
+    Q = s_Q[id]; dHxdr += D.x*Q; dHxds += D.y*Q; dHxdt += D.z*Q; id += BSIZE;
+    Q = s_Q[id]; dHydr += D.x*Q; dHyds += D.y*Q; dHydt += D.z*Q; id += BSIZE;
+    Q = s_Q[id]; dHzdr += D.x*Q; dHzds += D.y*Q; dHzdt += D.z*Q; id += BSIZE;
+    Q = s_Q[id]; dExdr += D.x*Q; dExds += D.y*Q; dExdt += D.z*Q; id += BSIZE;
+    Q = s_Q[id]; dEydr += D.x*Q; dEyds += D.y*Q; dEydt += D.z*Q; id += BSIZE;
+    Q = s_Q[id]; dEzdr += D.x*Q; dEzds += D.y*Q; dEzdt += D.z*Q; 
+
+    ++m;
+#if ( (p_Np) % 2 )==0
+    D = tex1Dfetch(t_DrDsDt, n+m*BSIZE);
+
+    id = m;
+    Q = s_Q[id]; dHxdr += D.x*Q; dHxds += D.y*Q; dHxdt += D.z*Q; id += BSIZE;
+    Q = s_Q[id]; dHydr += D.x*Q; dHyds += D.y*Q; dHydt += D.z*Q; id += BSIZE;
+    Q = s_Q[id]; dHzdr += D.x*Q; dHzds += D.y*Q; dHzdt += D.z*Q; id += BSIZE;
+    Q = s_Q[id]; dExdr += D.x*Q; dExds += D.y*Q; dExdt += D.z*Q; id += BSIZE;
+    Q = s_Q[id]; dEydr += D.x*Q; dEyds += D.y*Q; dEydt += D.z*Q; id += BSIZE;
+    Q = s_Q[id]; dEzdr += D.x*Q; dEzds += D.y*Q; dEzdt += D.z*Q; 
+
+    ++m;
+
+#if ( (p_Np)%3 )==0
+    D = tex1Dfetch(t_DrDsDt, n+m*BSIZE);
+
+    id = m;
+    Q = s_Q[id]; dHxdr += D.x*Q; dHxds += D.y*Q; dHxdt += D.z*Q; id += BSIZE;
+    Q = s_Q[id]; dHydr += D.x*Q; dHyds += D.y*Q; dHydt += D.z*Q; id += BSIZE;
+    Q = s_Q[id]; dHzdr += D.x*Q; dHzds += D.y*Q; dHzdt += D.z*Q; id += BSIZE;
+    Q = s_Q[id]; dExdr += D.x*Q; dExds += D.y*Q; dExdt += D.z*Q; id += BSIZE;
+    Q = s_Q[id]; dEydr += D.x*Q; dEyds += D.y*Q; dEydt += D.z*Q; id += BSIZE;
+    Q = s_Q[id]; dEzdr += D.x*Q; dEzds += D.y*Q; dEzdt += D.z*Q; 
+
+    ++m;
+#endif
+#endif
+  }
+  
+  const float drdx= s_facs[0];
+  const float drdy= s_facs[1];
+  const float drdz= s_facs[2];
+  const float dsdx= s_facs[4];
+  const float dsdy= s_facs[5];
+  const float dsdz= s_facs[6];
+  const float dtdx= s_facs[8];
+  const float dtdy= s_facs[9];
+  const float dtdz= s_facs[10];
+  
+  m = n+p_Nfields*BSIZE*k;
+
+  g_rhsQ[m] = -(drdy*dEzdr+dsdy*dEzds+dtdy*dEzdt - drdz*dEydr-dsdz*dEyds-dtdz*dEydt); m += BSIZE;
+  g_rhsQ[m] = -(drdz*dExdr+dsdz*dExds+dtdz*dExdt - drdx*dEzdr-dsdx*dEzds-dtdx*dEzdt); m += BSIZE;
+  g_rhsQ[m] = -(drdx*dEydr+dsdx*dEyds+dtdx*dEydt - drdy*dExdr-dsdy*dExds-dtdy*dExdt); m += BSIZE;
+  g_rhsQ[m] =  (drdy*dHzdr+dsdy*dHzds+dtdy*dHzdt - drdz*dHydr-dsdz*dHyds-dtdz*dHydt); m += BSIZE;
+  g_rhsQ[m] =  (drdz*dHxdr+dsdz*dHxds+dtdz*dHxdt - drdx*dHzdr-dsdx*dHzds-dtdx*dHzdt); m += BSIZE;
+  g_rhsQ[m] =  (drdx*dHydr+dsdx*dHyds+dtdx*dHydt - drdy*dHxdr-dsdy*dHxds-dtdy*dHxdt); 
+}
+
+__global__ void MaxwellsGPU_SURF_Kernel3D(float *g_Q, float *g_rhsQ){
+
+  __device__ __shared__ float s_fluxQ[p_Nfields*p_Nfp*p_Nfaces];
+
+  const int n = threadIdx.x;
+  const int k = blockIdx.x;
+  int m;
+
+  /* grab surface nodes and store flux in shared memory */
+  if(n< (p_Nfp*p_Nfaces) ){
+    /* coalesced reads (maybe) */
+    m = 7*(k*p_Nfp*p_Nfaces)+n;
+    const  int idM   = tex1Dfetch(t_surfinfo, m); m += p_Nfp*p_Nfaces;
+           int idP   = tex1Dfetch(t_surfinfo, m); m += p_Nfp*p_Nfaces;
+    const  float Fsc = tex1Dfetch(t_surfinfo, m); m += p_Nfp*p_Nfaces;
+    const  float Bsc = tex1Dfetch(t_surfinfo, m); m += p_Nfp*p_Nfaces;
+    const  float nx  = tex1Dfetch(t_surfinfo, m); m += p_Nfp*p_Nfaces;
+    const  float ny  = tex1Dfetch(t_surfinfo, m); m += p_Nfp*p_Nfaces;
+    const  float nz  = tex1Dfetch(t_surfinfo, m);
+
+    /* check if idP<0  */
+    double dHx, dHy, dHz, dEx, dEy, dEz;
+    if(idP<0){
+      idP = p_Nfields*(-1-idP);
+      
+      dHx = Fsc*(tex1Dfetch(t_partQ, idP+0) - tex1Dfetch(t_Q, idM+0*BSIZE));
+      dHy = Fsc*(tex1Dfetch(t_partQ, idP+1) - tex1Dfetch(t_Q, idM+1*BSIZE));
+      dHz = Fsc*(tex1Dfetch(t_partQ, idP+2) - tex1Dfetch(t_Q, idM+2*BSIZE));
+      
+      dEx = Fsc*(tex1Dfetch(t_partQ, idP+3) - tex1Dfetch(t_Q, idM+3*BSIZE));
+      dEy = Fsc*(tex1Dfetch(t_partQ, idP+4) - tex1Dfetch(t_Q, idM+4*BSIZE));
+      dEz = Fsc*(tex1Dfetch(t_partQ, idP+5) - tex1Dfetch(t_Q, idM+5*BSIZE));
+    }
+    else{
+      dHx = Fsc*(tex1Dfetch(t_Q, idP+0*BSIZE) - tex1Dfetch(t_Q, idM+0*BSIZE));
+      dHy = Fsc*(tex1Dfetch(t_Q, idP+1*BSIZE) - tex1Dfetch(t_Q, idM+1*BSIZE));
+      dHz = Fsc*(tex1Dfetch(t_Q, idP+2*BSIZE) - tex1Dfetch(t_Q, idM+2*BSIZE));
+      
+      dEx = Fsc*(Bsc*tex1Dfetch(t_Q, idP+3*BSIZE) - tex1Dfetch(t_Q, idM+3*BSIZE));
+      dEy = Fsc*(Bsc*tex1Dfetch(t_Q, idP+4*BSIZE) - tex1Dfetch(t_Q, idM+4*BSIZE));
+      dEz = Fsc*(Bsc*tex1Dfetch(t_Q, idP+5*BSIZE) - tex1Dfetch(t_Q, idM+5*BSIZE));
+    }
+
+    const double ndotdH = nx*dHx + ny*dHy + nz*dHz;
+    const double ndotdE = nx*dEx + ny*dEy + nz*dEz;
+
+    m = n;
+    s_fluxQ[m] = -ny*dEz + nz*dEy + dHx - ndotdH*nx; m += p_Nfp*p_Nfaces;
+    s_fluxQ[m] = -nz*dEx + nx*dEz + dHy - ndotdH*ny; m += p_Nfp*p_Nfaces;
+    s_fluxQ[m] = -nx*dEy + ny*dEx + dHz - ndotdH*nz; m += p_Nfp*p_Nfaces;
+
+    s_fluxQ[m] =  ny*dHz - nz*dHy + dEx - ndotdE*nx; m += p_Nfp*p_Nfaces;
+    s_fluxQ[m] =  nz*dHx - nx*dHz + dEy - ndotdE*ny; m += p_Nfp*p_Nfaces;
+    s_fluxQ[m] =  nx*dHy - ny*dHx + dEz - ndotdE*nz; 
+  }
+
+  /* make sure all element data points are cached */
+  __syncthreads();
+
+  if(n< (p_Np))
+  {
+    float rhsHx = 0, rhsHy = 0, rhsHz = 0;
+    float rhsEx = 0, rhsEy = 0, rhsEz = 0;
+    
+    int sk = n;
+    /* can manually unroll to 4 because there are 4 faces */
+    for(m=0;p_Nfaces*p_Nfp-m;){
+      const float4 L = tex1Dfetch(t_LIFT, sk); sk+=p_Np;
+
+      /* broadcast */
+      int sk1 = m;
+      rhsHx += L.x*s_fluxQ[sk1]; sk1 += p_Nfp*p_Nfaces;
+      rhsHy += L.x*s_fluxQ[sk1]; sk1 += p_Nfp*p_Nfaces;
+      rhsHz += L.x*s_fluxQ[sk1]; sk1 += p_Nfp*p_Nfaces;
+      rhsEx += L.x*s_fluxQ[sk1]; sk1 += p_Nfp*p_Nfaces;
+      rhsEy += L.x*s_fluxQ[sk1]; sk1 += p_Nfp*p_Nfaces;
+      rhsEz += L.x*s_fluxQ[sk1]; sk1 += p_Nfp*p_Nfaces;
+      ++m;
+
+      /* broadcast */
+      sk1 = m;
+      rhsHx += L.y*s_fluxQ[sk1]; sk1 += p_Nfp*p_Nfaces;
+      rhsHy += L.y*s_fluxQ[sk1]; sk1 += p_Nfp*p_Nfaces;
+      rhsHz += L.y*s_fluxQ[sk1]; sk1 += p_Nfp*p_Nfaces;
+      rhsEx += L.y*s_fluxQ[sk1]; sk1 += p_Nfp*p_Nfaces;
+      rhsEy += L.y*s_fluxQ[sk1]; sk1 += p_Nfp*p_Nfaces;
+      rhsEz += L.y*s_fluxQ[sk1]; sk1 += p_Nfp*p_Nfaces;
+      ++m;
+
+      /* broadcast */
+      sk1 = m;
+      rhsHx += L.z*s_fluxQ[sk1]; sk1 += p_Nfp*p_Nfaces;
+      rhsHy += L.z*s_fluxQ[sk1]; sk1 += p_Nfp*p_Nfaces;
+      rhsHz += L.z*s_fluxQ[sk1]; sk1 += p_Nfp*p_Nfaces;
+      rhsEx += L.z*s_fluxQ[sk1]; sk1 += p_Nfp*p_Nfaces;
+      rhsEy += L.z*s_fluxQ[sk1]; sk1 += p_Nfp*p_Nfaces;
+      rhsEz += L.z*s_fluxQ[sk1]; sk1 += p_Nfp*p_Nfaces;
+      ++m;
+
+      /* broadcast */
+      sk1 = m;
+      rhsHx += L.w*s_fluxQ[sk1]; sk1 += p_Nfp*p_Nfaces;
+      rhsHy += L.w*s_fluxQ[sk1]; sk1 += p_Nfp*p_Nfaces;
+      rhsHz += L.w*s_fluxQ[sk1]; sk1 += p_Nfp*p_Nfaces;
+      rhsEx += L.w*s_fluxQ[sk1]; sk1 += p_Nfp*p_Nfaces;
+      rhsEy += L.w*s_fluxQ[sk1]; sk1 += p_Nfp*p_Nfaces;
+      rhsEz += L.w*s_fluxQ[sk1]; sk1 += p_Nfp*p_Nfaces;
+      ++m;
+
+    }
+    
+    m = n+p_Nfields*k*BSIZE;
+    g_rhsQ[m] += rhsHx; m += BSIZE;
+    g_rhsQ[m] += rhsHy; m += BSIZE;
+    g_rhsQ[m] += rhsHz; m += BSIZE;
+    g_rhsQ[m] += rhsEx; m += BSIZE;
+    g_rhsQ[m] += rhsEy; m += BSIZE;
+    g_rhsQ[m] += rhsEz; m += BSIZE;
+
+  }
+}
+
+
+__global__ void MaxwellsGPU_RK_Kernel3D(int Ntotal, float *g_resQ, float *g_rhsQ, float *g_Q, float fa, float fb, float fdt){
+  
+  int n = blockIdx.x * blockDim.x + threadIdx.x;
+    
+  if(n<Ntotal){
+    float rhs = g_rhsQ[n];
+    float res = g_resQ[n];
+    res = fa*res + fdt*rhs;
+    
+    g_resQ[n] = res;
+    g_Q[n]    += fb*res;
+  }
+
+} 
+
+
+/* assumes data resides on device */
+void MaxwellsKernel3d(Mesh *mesh, float frka, float frkb, float fdt){
+
+  /* grab data from device and initiate sends */
+  MaxwellsMPISend3d(mesh);
+   
+  int ThreadsPerBlock, BlocksPerGrid;	
+
+  BlocksPerGrid    = mesh->K; 
+  ThreadsPerBlock = p_Np;
+
+  /* evaluate volume derivatives */
+  MaxwellsGPU_VOL_Kernel3D <<< BlocksPerGrid, ThreadsPerBlock >>>  (c_rhsQ);
+
+  /* finalize sends and recvs, and transfer to device */
+  MaxwellsMPIRecv3d(mesh, c_partQ);
+
+  BlocksPerGrid = mesh->K;
+
+  if( ( p_Nfp*p_Nfaces ) > (p_Np) )
+    ThreadsPerBlock = p_Nfp*p_Nfaces;
+  else
+    ThreadsPerBlock = p_Np;
+
+  /* evaluate surface contributions */
+  MaxwellsGPU_SURF_Kernel3D <<< BlocksPerGrid, ThreadsPerBlock >>> (c_Q, c_rhsQ);
+
+  int Ntotal = mesh->K*BSIZE*p_Nfields;
+  
+  ThreadsPerBlock = 256;
+  BlocksPerGrid = (Ntotal+ThreadsPerBlock-1)/ThreadsPerBlock;
+
+  /* update RK Step */
+  MaxwellsGPU_RK_Kernel3D<<< BlocksPerGrid, ThreadsPerBlock >>> 
+	  (Ntotal, c_resQ, c_rhsQ, c_Q, frka, frkb, fdt);
+
+}
+
+
+
+
+void gpu_set_data3d(int K,
+		  double *d_Hx, double *d_Hy, double *d_Hz,
+		  double *d_Ex, double *d_Ey, double *d_Ez){
+
+
+  float *f_Q = (float*) calloc(K*p_Nfields*BSIZE,sizeof(float));
+  
+  /* also load into usual data matrices */
+  
+  for(int k=0;k<K;++k){
+    int gk = k;
+    for(int n=0;n<p_Np;++n)
+      f_Q[n        +k*BSIZE*p_Nfields] = d_Hx[n+gk*p_Np];
+    for(int n=0;n<p_Np;++n)
+      f_Q[n  +BSIZE+k*BSIZE*p_Nfields] = d_Hy[n+gk*p_Np];
+    for(int n=0;n<p_Np;++n)
+      f_Q[n+2*BSIZE+k*BSIZE*p_Nfields] = d_Hz[n+gk*p_Np];
+    for(int n=0;n<p_Np;++n)
+      f_Q[n+3*BSIZE+k*BSIZE*p_Nfields] = d_Ex[n+gk*p_Np];
+    for(int n=0;n<p_Np;++n)
+      f_Q[n+4*BSIZE+k*BSIZE*p_Nfields] = d_Ey[n+gk*p_Np];
+    for(int n=0;n<p_Np;++n)
+      f_Q[n+5*BSIZE+k*BSIZE*p_Nfields] = d_Ez[n+gk*p_Np];
+  }
+  
+  cudaMemcpy(c_Q, f_Q, BSIZE*K*p_Nfields*sizeof(float), cudaMemcpyHostToDevice);
+  
+  free(f_Q);
+}
+  
+void gpu_get_data3d(int K,
+		  double *d_Hx, double *d_Hy, double *d_Hz,
+		  double *d_Ex, double *d_Ey, double *d_Ez){
+
+  float *f_Q = (float*) calloc(K*p_Nfields*BSIZE,sizeof(float));
+  
+  cudaMemcpy(f_Q, c_Q, K*BSIZE*p_Nfields*sizeof(float), cudaMemcpyDeviceToHost);
+
+  /* also load into usual data matrices */
+  
+  for(int k=0;k<K;++k){
+    int gk = k;
+    for(int n=0;n<p_Np;++n)
+      d_Hx[n+gk*p_Np] = f_Q[n        +k*BSIZE*p_Nfields];
+    for(int n=0;n<p_Np;++n) 
+      d_Hy[n+gk*p_Np] = f_Q[n  +BSIZE+k*BSIZE*p_Nfields];
+    for(int n=0;n<p_Np;++n)
+      d_Hz[n+gk*p_Np] = f_Q[n+2*BSIZE+k*BSIZE*p_Nfields];
+    for(int n=0;n<p_Np;++n)
+      d_Ex[n+gk*p_Np] = f_Q[n+3*BSIZE+k*BSIZE*p_Nfields];
+    for(int n=0;n<p_Np;++n) 
+      d_Ey[n+gk*p_Np] = f_Q[n+4*BSIZE+k*BSIZE*p_Nfields];
+    for(int n=0;n<p_Np;++n)
+      d_Ez[n+gk*p_Np] = f_Q[n+5*BSIZE+k*BSIZE*p_Nfields];
+
+  }
+
+  free(f_Q);
+}
+
+__global__ void partial_get_kernel3d(int Ntotal, int *g_index, float *g_partQ){
+  
+  int n = blockIdx.x * blockDim.x + threadIdx.x;
+    
+  if(n<Ntotal)
+    g_partQ[n] = tex1Dfetch(t_Q, g_index[n]);
+  
+} 
+
+void get_partial_gpu_data3d(int Ntotal, int *g_index, float *h_partQ){
+
+  int ThreadsPerBlock = 256;
+  int BlocksPerGrid = (Ntotal+ThreadsPerBlock-1)/ThreadsPerBlock;
+
+  partial_get_kernel3d <<< BlocksPerGrid, ThreadsPerBlock >>> (Ntotal, g_index, c_tmp);
+
+  cudaMemcpy(h_partQ, c_tmp, Ntotal*sizeof(float), cudaMemcpyDeviceToHost);
+}
+
+}