integrating recent changes from fermi-test into fermi

(i'll use "fermi" for more disruptive changes to the pipeline model such
as updating the MSHRs and getting rid of the warp tracker, ripping out DWF, etc...)

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

1 files changed, 0 insertions, 551 deletions

diff --git a/benchmarks/CUDA/DG/src/MaxwellsKernel3d.cu b/benchmarks/CUDA/DG/src/MaxwellsKernel3d.cu
deleted file mode 100644
index bfa4396..0000000
--- a/benchmarks/CUDA/DG/src/MaxwellsKernel3d.cu
+++ /dev/null
@@ -1,551 +0,0 @@
-/* -*- 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);
-}
-
-}