From d9eb72d78ce75110fb3439d68c0ed35808382e87 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Andr=C3=A9=20R=2E=20Brodtkorb?= <Andre.Brodtkorb@sintef.no>
Date: Thu, 1 Nov 2018 20:43:28 +0100
Subject: [PATCH] Refactoring

---
 GPUSimulators/Common.py                  |  11 +-
 GPUSimulators/EE2D_KP07_dimsplit.py      |  10 +-
 GPUSimulators/Simulator.py               |  13 +-
 GPUSimulators/cuda/EE2D_KP07_dimsplit.cu | 175 ++++++++++-----------
 GPUSimulators/cuda/common.h              | 192 +++++++++++------------
 GPUSimulators/cuda/limiters.h            |  40 ++---
 6 files changed, 220 insertions(+), 221 deletions(-)

diff --git a/GPUSimulators/Common.py b/GPUSimulators/Common.py
index b973692..a0c2fb3 100644
--- a/GPUSimulators/Common.py
+++ b/GPUSimulators/Common.py
@@ -92,7 +92,7 @@ class CudaArray2D:
         
         #self.logger.debug("Allocating [%dx%d] buffer", self.nx, self.ny)
         #Should perhaps use pycuda.driver.mem_alloc_data.pitch() here
-        self.data = pycuda.gpuarray.empty((ny_halo, nx_halo), dtype)
+        self.data = pycuda.gpuarray.zeros((ny_halo, nx_halo), dtype)
         
         #If we don't have any data, just allocate and return
         if cpu_data is None:
@@ -310,5 +310,12 @@ class ArakawaA2D:
         stream.synchronize()
         return cpu_variables
         
-        
+    """
+    Checks that data is still sane
+    """
+    def check(self):
+        for i, gpu_variable in enumerate(self.gpu_variables):
+            var_sum = pycuda.gpuarray.sum(gpu_variable.data).get()
+            self.logger.debug("Data %d with size [%d x %d] has sum %f", i, gpu_variable.nx, gpu_variable.ny, var_sum)
+            assert np.isnan(var_sum) == False, "Data contains NaN values!"
     
\ No newline at end of file
diff --git a/GPUSimulators/EE2D_KP07_dimsplit.py b/GPUSimulators/EE2D_KP07_dimsplit.py
index 84fcc33..99ba8c7 100644
--- a/GPUSimulators/EE2D_KP07_dimsplit.py
+++ b/GPUSimulators/EE2D_KP07_dimsplit.py
@@ -57,7 +57,7 @@ class EE2D_KP07_dimsplit (Simulator.BaseSimulator):
                  dx, dy, dt, \
                  gamma, \
                  theta=1.3, \
-                 block_width=8, block_height=4):
+                 block_width=16, block_height=8):
                  
         # Call super constructor
         super().__init__(context, \
@@ -65,10 +65,8 @@ class EE2D_KP07_dimsplit (Simulator.BaseSimulator):
             dx, dy, dt, \
             block_width, block_height)
         self.gamma = np.float32(gamma)
-            
         self.theta = np.float32(theta)
 
-        #Get kernels
         #Get kernels
         self.kernel = context.get_prepared_kernel("cuda/EE2D_KP07_dimsplit.cu", "KP07DimsplitKernel", \
                                         "iifffffiPiPiPiPiPiPiPiPi", \
@@ -135,4 +133,8 @@ class EE2D_KP07_dimsplit (Simulator.BaseSimulator):
         self.t += dt
         
     def download(self):
-        return self.u0.download(self.stream)
\ No newline at end of file
+        return self.u0.download(self.stream)
+
+    def check(self):
+        self.u0.check()
+        self.u1.check()
\ No newline at end of file
diff --git a/GPUSimulators/Simulator.py b/GPUSimulators/Simulator.py
index 3792a01..94ef35b 100644
--- a/GPUSimulators/Simulator.py
+++ b/GPUSimulators/Simulator.py
@@ -122,6 +122,7 @@ class BaseSimulator:
                 if (t.elapsed() >= next_print):
                     self.logger.info("%s simulated %d of %d steps (Euler)", self, i, n)
                     next_print += self.log_every
+                    self.check()
     
             
         self.logger.info("%s simulated %f seconds to %f with %d steps (Euler)", self, t_end, self.t, n)
@@ -153,6 +154,7 @@ class BaseSimulator:
                 if (t.elapsed() >= next_print):
                     self.logger.info("%s simulated %d of %d steps (RK2)", self, i, n)
                     next_print += self.log_every
+                    self.check()
             
         self.logger.info("%s simulated %f seconds to %f with %d steps (RK2)", self, t_end, self.t, n)
         return self.t, n
@@ -184,6 +186,7 @@ class BaseSimulator:
                 if (t.elapsed() >= next_print):
                     self.logger.info("%s simulated %d of %d steps (Dimsplit)", self, i, n)
                     next_print += self.log_every
+                    self.check()
             
         self.logger.info("%s simulated %f seconds to %f with %d steps (Dimsplit)", self, t_end, self.t, 2*n)
         return self.t, 2*n
@@ -203,13 +206,17 @@ class BaseSimulator:
         
     def stepDimsplitYX(self, dt):
         raise(NotImplementedError("Needs to be implemented in subclass"))
-        
-    def sim_time(self):
-        return self.t
 
     def download(self):
         raise(NotImplementedError("Needs to be implemented in subclass"))
+
+    def check(self):
+        raise(NotImplementedError("Needs to be implemented in subclass"))
+        
+    def sim_time(self):
+        return self.t
         
     def synchronize(self):
         self.stream.synchronize()
 
+        
diff --git a/GPUSimulators/cuda/EE2D_KP07_dimsplit.cu b/GPUSimulators/cuda/EE2D_KP07_dimsplit.cu
index 370efbd..2397712 100644
--- a/GPUSimulators/cuda/EE2D_KP07_dimsplit.cu
+++ b/GPUSimulators/cuda/EE2D_KP07_dimsplit.cu
@@ -27,97 +27,95 @@ along with this program.  If not, see <http://www.gnu.org/licenses/>.
 
 __device__
 void computeFluxF(float Q[4][BLOCK_HEIGHT+4][BLOCK_WIDTH+4],
-                  float Qx[4][BLOCK_HEIGHT+2][BLOCK_WIDTH+2],
-                  float F[4][BLOCK_HEIGHT+1][BLOCK_WIDTH+1],
-                  const float gamma_, const float dx_, const float dt_) {    
-    int j=threadIdx.y;
-    const int l = j + 2; //Skip ghost cells
-    for (int i=threadIdx.x; i<BLOCK_WIDTH+1; i+=BLOCK_WIDTH) {
-        const int k = i + 1;
-        // Reconstruct point values of Q at the left and right hand side 
-        // of the cell for both the left (i) and right (i+1) cell 
-        const float4 Q_rl = make_float4(Q[0][l][k+1] - 0.5f*Qx[0][j][i+1],
-                                        Q[1][l][k+1] - 0.5f*Qx[1][j][i+1],
-                                        Q[2][l][k+1] - 0.5f*Qx[2][j][i+1],
-                                        Q[3][l][k+1] - 0.5f*Qx[3][j][i+1]);
-        const float4 Q_rr = make_float4(Q[0][l][k+1] + 0.5f*Qx[0][j][i+1],
-                                        Q[1][l][k+1] + 0.5f*Qx[1][j][i+1],
-                                        Q[2][l][k+1] + 0.5f*Qx[2][j][i+1],
-                                        Q[3][l][k+1] + 0.5f*Qx[3][j][i+1]);
-                                     
-        const float4 Q_ll = make_float4(Q[0][l][k] - 0.5f*Qx[0][j][i],
-                                        Q[1][l][k] - 0.5f*Qx[1][j][i],
-                                        Q[2][l][k] - 0.5f*Qx[2][j][i],
-                                        Q[3][l][k] - 0.5f*Qx[3][j][i]);
-        const float4 Q_lr = make_float4(Q[0][l][k] + 0.5f*Qx[0][j][i],
-                                        Q[1][l][k] + 0.5f*Qx[1][j][i],
-                                        Q[2][l][k] + 0.5f*Qx[2][j][i],
-                                        Q[3][l][k] + 0.5f*Qx[3][j][i]);
-                                
-        //Evolve half a timestep (predictor step)
-        const float4 Q_r_bar = Q_rl + dt_/(2.0f*dx_) * (F_func(Q_rl, gamma_) - F_func(Q_rr, gamma_));
-        const float4 Q_l_bar = Q_lr + dt_/(2.0f*dx_) * (F_func(Q_ll, gamma_) - F_func(Q_lr, gamma_));
+                  float Qx[4][BLOCK_HEIGHT+4][BLOCK_WIDTH+4],
+                  float F[4][BLOCK_HEIGHT+4][BLOCK_WIDTH+4],
+                  const float gamma_, const float dx_, const float dt_) {
+    for (int j=threadIdx.y; j<BLOCK_HEIGHT+4; j+=BLOCK_HEIGHT) {
+        for (int i=threadIdx.x+1; i<BLOCK_WIDTH+2; i+=BLOCK_WIDTH) {
+            // Reconstruct point values of Q at the left and right hand side 
+            // of the cell for both the left (i) and right (i+1) cell 
+            const float4 Q_rl = make_float4(Q[0][j][i+1] - 0.5f*Qx[0][j][i+1],
+                                            Q[1][j][i+1] - 0.5f*Qx[1][j][i+1],
+                                            Q[2][j][i+1] - 0.5f*Qx[2][j][i+1],
+                                            Q[3][j][i+1] - 0.5f*Qx[3][j][i+1]);
+            const float4 Q_rr = make_float4(Q[0][j][i+1] + 0.5f*Qx[0][j][i+1],
+                                            Q[1][j][i+1] + 0.5f*Qx[1][j][i+1],
+                                            Q[2][j][i+1] + 0.5f*Qx[2][j][i+1],
+                                            Q[3][j][i+1] + 0.5f*Qx[3][j][i+1]);
 
-        // Compute flux based on prediction
-        const float4 flux = CentralUpwindFlux(Q_l_bar, Q_r_bar, gamma_);
-        
-        //Write to shared memory
-        F[0][j][i] = flux.x;
-        F[1][j][i] = flux.y;
-        F[2][j][i] = flux.z;
-        F[3][j][i] = flux.w;
+            const float4 Q_ll = make_float4(Q[0][j][i] - 0.5f*Qx[0][j][i],
+                                            Q[1][j][i] - 0.5f*Qx[1][j][i],
+                                            Q[2][j][i] - 0.5f*Qx[2][j][i],
+                                            Q[3][j][i] - 0.5f*Qx[3][j][i]);
+            const float4 Q_lr = make_float4(Q[0][j][i] + 0.5f*Qx[0][j][i],
+                                            Q[1][j][i] + 0.5f*Qx[1][j][i],
+                                            Q[2][j][i] + 0.5f*Qx[2][j][i],
+                                            Q[3][j][i] + 0.5f*Qx[3][j][i]);
+
+
+            //Evolve half a timestep (predictor step)
+            const float4 Q_r_bar = Q_rl + dt_/(2.0f*dx_) * (F_func(Q_rl, gamma_) - F_func(Q_rr, gamma_));
+            const float4 Q_l_bar = Q_lr + dt_/(2.0f*dx_) * (F_func(Q_ll, gamma_) - F_func(Q_lr, gamma_));
+
+            // Compute flux based on prediction
+            const float4 flux = CentralUpwindFlux(Q_l_bar, Q_r_bar, gamma_);
+            
+            //Write to shared memory
+            F[0][j][i] = flux.x;
+            F[1][j][i] = flux.y;
+            F[2][j][i] = flux.z;
+            F[3][j][i] = flux.w;
+        }
     }
 }
 
 __device__
 void computeFluxG(float Q[4][BLOCK_HEIGHT+4][BLOCK_WIDTH+4],
-                  float Qy[4][BLOCK_HEIGHT+2][BLOCK_WIDTH+2],
-                  float G[4][BLOCK_HEIGHT+1][BLOCK_WIDTH+1],
-                  const float gamma_, const float dy_, const float dt_) {    
-    int i=threadIdx.x;
-    const int k = i + 2; //Skip ghost cells
-    for (int j=threadIdx.y; j<BLOCK_HEIGHT+1; j+=BLOCK_HEIGHT) {
-        const int l = j + 1;
-        // Reconstruct point values of Q at the left and right hand side 
-        // of the cell for both the left (i) and right (i+1) cell 
-        //NOte that hu and hv are swapped ("transposing" the domain)!
-        const float4 Q_rl = make_float4(Q[0][l+1][k] - 0.5f*Qy[0][j+1][i],
-                                        Q[2][l+1][k] - 0.5f*Qy[2][j+1][i],
-                                        Q[1][l+1][k] - 0.5f*Qy[1][j+1][i],
-                                        Q[3][l+1][k] - 0.5f*Qy[3][j+1][i]);
-        const float4 Q_rr = make_float4(Q[0][l+1][k] + 0.5f*Qy[0][j+1][i],
-                                        Q[2][l+1][k] + 0.5f*Qy[2][j+1][i],
-                                        Q[1][l+1][k] + 0.5f*Qy[1][j+1][i],
-                                        Q[3][l+1][k] + 0.5f*Qy[3][j+1][i]);
+                  float Qy[4][BLOCK_HEIGHT+4][BLOCK_WIDTH+4],
+                  float G[4][BLOCK_HEIGHT+4][BLOCK_WIDTH+4],
+                  const float gamma_, const float dy_, const float dt_) {
+    for (int j=threadIdx.y+1; j<BLOCK_HEIGHT+2; j+=BLOCK_HEIGHT) {
+        for (int i=threadIdx.x; i<BLOCK_WIDTH+4; i+=BLOCK_WIDTH) {
+            // Reconstruct point values of Q at the left and right hand side 
+            // of the cell for both the left (i) and right (i+1) cell 
+            //NOte that hu and hv are swapped ("transposing" the domain)!
+            const float4 Q_rl = make_float4(Q[0][j+1][i] - 0.5f*Qy[0][j+1][i],
+                                            Q[2][j+1][i] - 0.5f*Qy[2][j+1][i],
+                                            Q[1][j+1][i] - 0.5f*Qy[1][j+1][i],
+                                            Q[3][j+1][i] - 0.5f*Qy[3][j+1][i]);
+            const float4 Q_rr = make_float4(Q[0][j+1][i] + 0.5f*Qy[0][j+1][i],
+                                            Q[2][j+1][i] + 0.5f*Qy[2][j+1][i],
+                                            Q[1][j+1][i] + 0.5f*Qy[1][j+1][i],
+                                            Q[3][j+1][i] + 0.5f*Qy[3][j+1][i]);
 
-        const float4 Q_ll = make_float4(Q[0][l][k] - 0.5f*Qy[0][j][i],
-                                        Q[2][l][k] - 0.5f*Qy[2][j][i],
-                                        Q[1][l][k] - 0.5f*Qy[1][j][i],
-                                        Q[3][l][k] - 0.5f*Qy[3][j][i]);
-        const float4 Q_lr = make_float4(Q[0][l][k] + 0.5f*Qy[0][j][i],
-                                        Q[2][l][k] + 0.5f*Qy[2][j][i],
-                                        Q[1][l][k] + 0.5f*Qy[1][j][i],
-                                        Q[3][l][k] + 0.5f*Qy[3][j][i]);
+            const float4 Q_ll = make_float4(Q[0][j][i] - 0.5f*Qy[0][j][i],
+                                            Q[2][j][i] - 0.5f*Qy[2][j][i],
+                                            Q[1][j][i] - 0.5f*Qy[1][j][i],
+                                            Q[3][j][i] - 0.5f*Qy[3][j][i]);
+            const float4 Q_lr = make_float4(Q[0][j][i] + 0.5f*Qy[0][j][i],
+                                            Q[2][j][i] + 0.5f*Qy[2][j][i],
+                                            Q[1][j][i] + 0.5f*Qy[1][j][i],
+                                            Q[3][j][i] + 0.5f*Qy[3][j][i]);
 
-        //Evolve half a timestep (predictor step)
-        const float4 Q_r_bar = Q_rl + dt_/(2.0f*dy_) * (F_func(Q_rl, gamma_) - F_func(Q_rr, gamma_));
-        const float4 Q_l_bar = Q_lr + dt_/(2.0f*dy_) * (F_func(Q_ll, gamma_) - F_func(Q_lr, gamma_));
-        
-        // Compute flux based on prediction
-        const float4 flux = CentralUpwindFlux(Q_l_bar, Q_r_bar, gamma_);
-        
-        //Write to shared memory
-        //Note that we here swap hu and hv back to the original
-        G[0][j][i] = flux.x;
-        G[1][j][i] = flux.z;
-        G[2][j][i] = flux.y;
-        G[3][j][i] = flux.w;
+            //Evolve half a timestep (predictor step)
+            const float4 Q_r_bar = Q_rl + dt_/(2.0f*dy_) * (F_func(Q_rl, gamma_) - F_func(Q_rr, gamma_));
+            const float4 Q_l_bar = Q_lr + dt_/(2.0f*dy_) * (F_func(Q_ll, gamma_) - F_func(Q_lr, gamma_));
+            
+            // Compute flux based on prediction
+            const float4 flux = CentralUpwindFlux(Q_l_bar, Q_r_bar, gamma_);
+            
+            //Write to shared memory
+            //Note that we here swap hu and hv back to the original
+            G[0][j][i] = flux.x;
+            G[1][j][i] = flux.z;
+            G[2][j][i] = flux.y;
+            G[3][j][i] = flux.w;
+        }
     }
 }
 
 
 
-
 /**
   * This unsplit kernel computes the 2D numerical scheme with a TVD RK2 time integration scheme
   */
@@ -150,8 +148,8 @@ __global__ void KP07DimsplitKernel(
         
     //Shared memory variables
     __shared__ float  Q[4][h+4][w+4];
-    __shared__ float Qx[4][h+2][w+2];
-    __shared__ float  F[4][h+1][w+1];
+    __shared__ float Qx[4][h+4][w+4];
+    __shared__ float  F[4][h+4][w+4];
     
     
     
@@ -162,14 +160,12 @@ __global__ void KP07DimsplitKernel(
     readBlock<w, h, gc>(    E0_ptr_,     E0_pitch_, Q[3], nx_+4, ny_+4);
     __syncthreads();
     
-    
     //Fix boundary conditions
     noFlowBoundary<w, h, gc,  1,  1>(Q[0], nx_, ny_);
-    noFlowBoundary<w, h, gc, -1,  1>(Q[1], nx_, ny_);
+    noFlowBoundary<w, h, gc,  1,  1>(Q[1], nx_, ny_);
     noFlowBoundary<w, h, gc,  1, -1>(Q[2], nx_, ny_);
     noFlowBoundary<w, h, gc,  1,  1>(Q[3], nx_, ny_);
     __syncthreads();
-    
 
 
     //Step 0 => evolve x first, then y
@@ -186,12 +182,11 @@ __global__ void KP07DimsplitKernel(
 
         //Set boundary conditions
         noFlowBoundary<w, h, gc,  1,  1>(Q[0], nx_, ny_);
-        noFlowBoundary<w, h, gc, -1,  1>(Q[1], nx_, ny_);
+        noFlowBoundary<w, h, gc,  1,  1>(Q[1], nx_, ny_);
         noFlowBoundary<w, h, gc,  1, -1>(Q[2], nx_, ny_);
         noFlowBoundary<w, h, gc,  1,  1>(Q[3], nx_, ny_);
         __syncthreads();
 
-
         //Compute fluxes along the y axis and evolve
         minmodSlopeY<w, h, gc, vars>(Q, Qx, theta_);
         __syncthreads();
@@ -200,7 +195,7 @@ __global__ void KP07DimsplitKernel(
         __syncthreads();
 
         evolveG<w, h, gc, vars>(Q, F, dy_, dt_);
-        __syncthreads();
+        __syncthreads();    
 
     }
     //Step 1 => evolve y first, then x
@@ -208,14 +203,16 @@ __global__ void KP07DimsplitKernel(
         //Compute fluxes along the y axis and evolve
         minmodSlopeY<w, h, gc, vars>(Q, Qx, theta_);
         __syncthreads();
+  
         computeFluxG(Q, Qx, F, gamma_, dy_, dt_);
         __syncthreads();
+  
         evolveG<w, h, gc, vars>(Q, F, dy_, dt_);
         __syncthreads();
-        
+  
         //Set boundary conditions
         noFlowBoundary<w, h, gc,  1,  1>(Q[0], nx_, ny_);
-        noFlowBoundary<w, h, gc, -1,  1>(Q[1], nx_, ny_);
+        noFlowBoundary<w, h, gc,  1,  1>(Q[1], nx_, ny_);
         noFlowBoundary<w, h, gc,  1, -1>(Q[2], nx_, ny_);
         noFlowBoundary<w, h, gc,  1,  1>(Q[3], nx_, ny_);
         __syncthreads();
@@ -223,8 +220,10 @@ __global__ void KP07DimsplitKernel(
         //Compute fluxes along the x axis and evolve
         minmodSlopeX<w, h, gc, vars>(Q, Qx, theta_);
         __syncthreads();
+
         computeFluxF(Q, Qx, F, gamma_, dx_, dt_);
         __syncthreads();
+
         evolveF<w, h, gc, vars>(Q, F, dx_, dt_);
         __syncthreads();
     }
diff --git a/GPUSimulators/cuda/common.h b/GPUSimulators/cuda/common.h
index 4a251e3..f936fce 100644
--- a/GPUSimulators/cuda/common.h
+++ b/GPUSimulators/cuda/common.h
@@ -151,86 +151,91 @@ inline __device__ void writeBlock(float* ptr_, int pitch_,
 
 template<int block_width, int block_height, int ghost_cells, int scale_east_west=1, int scale_north_south=1>
 __device__ void noFlowBoundary(float Q[block_height+2*ghost_cells][block_width+2*ghost_cells], const int nx_, const int ny_) {
-    const int ti = blockDim.x*blockIdx.x + threadIdx.x + ghost_cells;
-    const int tj = blockDim.y*blockIdx.y + threadIdx.y + ghost_cells;
     
-    const int i = threadIdx.x + ghost_cells;
-    const int j = threadIdx.y + ghost_cells;
-    
-    
-    // West boundary
-    if (ti == ghost_cells) {
-        Q[j][i-1] = scale_east_west*Q[j][i];
-    }
-    if (ghost_cells >= 2 && ti == ghost_cells + 1) {
-        Q[j][i-3] = scale_east_west*Q[j][i];
-    }
-    if (ghost_cells >= 3 && ti == ghost_cells + 2) {
-        Q[j][i-5] = scale_east_west*Q[j][i];
-    }
-    if (ghost_cells >= 4 && ti == ghost_cells + 3) {
-        Q[j][i-7] = scale_east_west*Q[j][i];
-    }
-    if (ghost_cells >= 5 && ti == ghost_cells + 4) {
-        Q[j][i-9] = scale_east_west*Q[j][i];
+    for (int j=threadIdx.y; j<block_height+2*ghost_cells; j+= block_height) {
+        const int i = threadIdx.x + ghost_cells;
+        const int ti = blockDim.x*blockIdx.x + i;
+        const int tj = blockDim.y*blockIdx.y + j;
+        
+        // West boundary
+        if (ti == ghost_cells) {
+            Q[j][i-1] = scale_east_west*Q[j][i];
+        }
+        if (ghost_cells >= 2 && ti == ghost_cells + 1) {
+            Q[j][i-3] = scale_east_west*Q[j][i];
+        }
+        if (ghost_cells >= 3 && ti == ghost_cells + 2) {
+            Q[j][i-5] = scale_east_west*Q[j][i];
+        }
+        if (ghost_cells >= 4 && ti == ghost_cells + 3) {
+            Q[j][i-7] = scale_east_west*Q[j][i];
+        }
+        if (ghost_cells >= 5 && ti == ghost_cells + 4) {
+            Q[j][i-9] = scale_east_west*Q[j][i];
+        }
+        
+        
+        
+        // East boundary
+        if (ti == nx_ + ghost_cells - 1) {
+            Q[j][i+1] = scale_east_west*Q[j][i];
+        }
+        if (ghost_cells >= 2 && ti == nx_ + ghost_cells - 2) {
+            Q[j][i+3] = scale_east_west*Q[j][i];
+        }
+        if (ghost_cells >= 3 && ti == nx_ + ghost_cells - 3) {
+            Q[j][i+5] = scale_east_west*Q[j][i];
+        }
+        if (ghost_cells >= 4 && ti == nx_ + ghost_cells - 4) {
+            Q[j][i+7] = scale_east_west*Q[j][i];
+        }
+        if (ghost_cells >= 5 && ti == nx_ + ghost_cells - 5) {
+            Q[j][i+9] = scale_east_west*Q[j][i];
+        }
     }
     
     
-    
-    // East boundary
-    if (ti == nx_ + ghost_cells - 1) {
-        Q[j][i+1] = scale_east_west*Q[j][i];
-    }
-    if (ghost_cells >= 2 && ti == nx_ + ghost_cells - 2) {
-        Q[j][i+3] = scale_east_west*Q[j][i];
-    }
-    if (ghost_cells >= 3 && ti == nx_ + ghost_cells - 3) {
-        Q[j][i+5] = scale_east_west*Q[j][i];
-    }
-    if (ghost_cells >= 3 && ti == nx_ + ghost_cells - 4) {
-        Q[j][i+7] = scale_east_west*Q[j][i];
-    }
-    if (ghost_cells >= 3 && ti == nx_ + ghost_cells - 5) {
-        Q[j][i+9] = scale_east_west*Q[j][i];
-    }
-    
-    
-    
-    
-    // South boundary
-    if (tj == ghost_cells) {
-        Q[j-1][i] = scale_north_south*Q[j][i];
-    }
-    if (ghost_cells >= 2 && tj == ghost_cells + 1) {
-        Q[j-3][i] = scale_north_south*Q[j][i];
-    }
-    if (ghost_cells >= 3 && tj == ghost_cells + 2) {
-        Q[j-5][i] = scale_north_south*Q[j][i];
-    }
-    if (ghost_cells >= 4 && tj == ghost_cells + 3) {
-        Q[j-7][i] = scale_north_south*Q[j][i];
-    }
-    if (ghost_cells >= 5 && tj == ghost_cells + 4) {
-        Q[j-9][i] = scale_north_south*Q[j][i];
-    }
-    
-    
-    
-    // North boundary
-    if (tj == ny_ + ghost_cells - 1) {
-        Q[j+1][i] = scale_north_south*Q[j][i];
-    }
-    if (ghost_cells >= 2 && tj == ny_ + ghost_cells - 2) {
-        Q[j+3][i] = scale_north_south*Q[j][i];
-    }
-    if (ghost_cells >= 3 && tj == ny_ + ghost_cells - 3) {
-        Q[j+5][i] = scale_north_south*Q[j][i];
-    }
-    if (ghost_cells >= 3 && tj == ny_ + ghost_cells - 4) {
-        Q[j+7][i] = scale_north_south*Q[j][i];
-    }
-    if (ghost_cells >= 3 && tj == ny_ + ghost_cells - 5) {
-        Q[j+9][i] = scale_north_south*Q[j][i];
+
+    for (int i=threadIdx.x; i<block_width+2*ghost_cells; i+= block_width) {
+        const int j = threadIdx.y + ghost_cells;
+        const int ti = blockDim.x*blockIdx.x + i;
+        const int tj = blockDim.y*blockIdx.y + j;
+
+        // South boundary
+        if (tj == ghost_cells) {
+            Q[j-1][i] = scale_north_south*Q[j][i];
+        }
+        if (ghost_cells >= 2 && tj == ghost_cells + 1) {
+            Q[j-3][i] = scale_north_south*Q[j][i];
+        }
+        if (ghost_cells >= 3 && tj == ghost_cells + 2) {
+            Q[j-5][i] = scale_north_south*Q[j][i];
+        }
+        if (ghost_cells >= 4 && tj == ghost_cells + 3) {
+            Q[j-7][i] = scale_north_south*Q[j][i];
+        }
+        if (ghost_cells >= 5 && tj == ghost_cells + 4) {
+            Q[j-9][i] = scale_north_south*Q[j][i];
+        }
+        
+        
+        
+        // North boundary
+        if (tj == ny_ + ghost_cells - 1) {
+            Q[j+1][i] = scale_north_south*Q[j][i];
+        }
+        if (ghost_cells >= 2 && tj == ny_ + ghost_cells - 2) {
+            Q[j+3][i] = scale_north_south*Q[j][i];
+        }
+        if (ghost_cells >= 3 && tj == ny_ + ghost_cells - 3) {
+            Q[j+5][i] = scale_north_south*Q[j][i];
+        }
+        if (ghost_cells >= 4 && tj == ny_ + ghost_cells - 4) {
+            Q[j+7][i] = scale_north_south*Q[j][i];
+        }
+        if (ghost_cells >= 5 && tj == ny_ + ghost_cells - 5) {
+            Q[j+9][i] = scale_north_south*Q[j][i];
+        }
     }
 }
 
@@ -247,21 +252,14 @@ __device__ void noFlowBoundary(float Q[block_height+2*ghost_cells][block_width+2
 
 template<int block_width, int block_height, int ghost_cells, int vars>
 __device__ void evolveF(float Q[vars][block_height+2*ghost_cells][block_width+2*ghost_cells],
-              float F[vars][block_height+1][block_width+1],
+              float F[vars][block_height+2*ghost_cells][block_width+2*ghost_cells],
               const float dx_, const float dt_) {
-    //Index of thread within block
-    const int tx = threadIdx.x;
-    const int ty = threadIdx.y;
-        
-    const int i = tx + ghost_cells; //Skip local ghost cells
-    const int j = ty + ghost_cells;
-    
-    //Index of cell within domain
-    //const int ti = blockDim.x*blockIdx.x + threadIdx.x + ghost_cells; //Skip global ghost cells, i.e., +1
-    //const int tj = blockDim.y*blockIdx.y + threadIdx.y + ghost_cells;
-    //if (ti > ghost_cells-1 && ti < nx_+ghost_cells && tj > ghost_cells-1 && tj < ny_+ghost_cells) {
     for (int var=0; var < vars; ++var) {
-        Q[var][j][i] = Q[var][j][i] + (F[var][ty][tx] - F[var][ty][tx+1]) * dt_ / dx_;
+        for (int j=threadIdx.y; j<block_height+2*ghost_cells; j+=block_height) {
+            for (int i=threadIdx.x+1; i<block_width+2*ghost_cells; i+=block_width) {
+                Q[var][j][i] = Q[var][j][i] + (F[var][j][i-1] - F[var][j][i]) * dt_ / dx_;
+            }
+        }
     }
 }
 
@@ -275,17 +273,14 @@ __device__ void evolveF(float Q[vars][block_height+2*ghost_cells][block_width+2*
   */
 template<int block_width, int block_height, int ghost_cells, int vars>
 __device__ void evolveG(float Q[vars][block_height+2*ghost_cells][block_width+2*ghost_cells],
-              float G[vars][block_height+1][block_width+1],
+              float G[vars][block_height+2*ghost_cells][block_width+2*ghost_cells],
               const float dy_, const float dt_) {
-    //Index of thread within block
-    const int tx = threadIdx.x;
-    const int ty = threadIdx.y;
-    
-    const int i = tx + ghost_cells; //Skip local ghost cells, i.e., +1
-    const int j = ty + ghost_cells;
-    
     for (int var=0; var < vars; ++var) {
-        Q[var][j][i] = Q[var][j][i] + (G[var][ty][tx] - G[var][ty+1][tx]) * dt_ / dy_;
+        for (int j=threadIdx.y+1; j<block_height+2*ghost_cells; j+=block_height) {
+            for (int i=threadIdx.x; i<block_width+2*ghost_cells; i+=block_width) {
+                Q[var][j][i] = Q[var][j][i] + (G[var][j-1][i] - G[var][j][i]) * dt_ / dy_;
+            }
+        }
     }
 }
 
@@ -328,5 +323,6 @@ __device__ void memset(float Q[vars][shmem_height][shmem_width], float value) {
 
 
 
+
 
 
diff --git a/GPUSimulators/cuda/limiters.h b/GPUSimulators/cuda/limiters.h
index 8d438ef..65ed892 100644
--- a/GPUSimulators/cuda/limiters.h
+++ b/GPUSimulators/cuda/limiters.h
@@ -47,21 +47,15 @@ __device__ __inline__ float minmodSlope(float left, float center, float right, f
   * Reconstructs a minmod slope for a whole block along the abscissa
   */
 template<int block_width, int block_height, int ghost_cells, int vars>
-__device__ void minmodSlopeX(float  Q[vars][block_height+2*ghost_cells][block_width+2*ghost_cells],
-                  float Qx[vars][block_height+2*(ghost_cells-1)][block_width+2*(ghost_cells-1)],
+__device__ void minmodSlopeX(float Q[vars][block_height+2*ghost_cells][block_width+2*ghost_cells],
+                  float Qx[vars][block_height+2*ghost_cells][block_width+2*ghost_cells],
                   const float theta_) {
-    //Index of thread within block
-    const int tx = threadIdx.x;
-    const int ty = threadIdx.y;
-    
-    const int j = ty;
-    const int l = j + ghost_cells; //Skip ghost cells
-
     //Reconstruct slopes along x axis
-    for (int i=tx; i<block_width+2*(ghost_cells-1); i+=block_width) {
-        const int k = i + 1;
-        for (int p=0; p<vars; ++p) {
-            Qx[p][j][i] = minmodSlope(Q[p][l][k-1], Q[p][l][k], Q[p][l][k+1], theta_);
+    for (int p=0; p<vars; ++p) {
+        for (int j=threadIdx.y; j<block_height+2*ghost_cells; j+=block_height) {
+            for (int i=threadIdx.x+1; i<block_width+3; i+=block_width) {
+                Qx[p][j][i] = minmodSlope(Q[p][j][i-1], Q[p][j][i], Q[p][j][i+1], theta_);
+            }
         }
     }
 }
@@ -71,21 +65,15 @@ __device__ void minmodSlopeX(float  Q[vars][block_height+2*ghost_cells][block_wi
   * Reconstructs a minmod slope for a whole block along the ordinate
   */
 template<int block_width, int block_height, int ghost_cells, int vars>
-__device__ void minmodSlopeY(float  Q[vars][block_height+2*ghost_cells][block_width+2*ghost_cells],
-                  float Qy[vars][block_height+2*(ghost_cells-1)][block_width+2*(ghost_cells-1)],
+__device__ void minmodSlopeY(float Q[vars][block_height+2*ghost_cells][block_width+2*ghost_cells],
+                  float Qy[vars][block_height+2*ghost_cells][block_width+2*ghost_cells],
                   const float theta_) {
-    //Index of thread within block
-    const int tx = threadIdx.x;
-    const int ty = threadIdx.y;
-    
-    const int i = tx;
-    const int k = i + ghost_cells; //Skip ghost cells
-
     //Reconstruct slopes along y axis
-    for (int j=ty; j<block_height+2*(ghost_cells-1); j+=block_height) {
-        const int l = j + 1;
-        for (int p=0; p<vars; ++p) {
-            Qy[p][j][i] = minmodSlope(Q[p][l-1][k], Q[p][l][k], Q[p][l+1][k], theta_);
+    for (int p=0; p<vars; ++p) {
+        for (int j=threadIdx.y+1; j<block_height+3; j+=block_height) {
+            for (int i=threadIdx.x; i<block_width+2*ghost_cells; i+=block_width) {
+                Qy[p][j][i] = minmodSlope(Q[p][j-1][i], Q[p][j][i], Q[p][j+1][i], theta_);
+            }
         }
     }
 }