 #include <stdlib.h>
 #include <stdio.h>
 #include <math.h>
-#include <time.h>
+
+ #include "timer.h"
 
 // CUDA libraries
 #include <cuda.h>
-#include <cutil.h>
+#include <helper_cuda.h>
 
 /** Problem size along one side; total number of cells is this squared */
 #define NUM 256
 
 // block size
-#define BLOCK_SIZE 64
+#define BLOCK_SIZE 128
 
 /** Double precision */
-//#define DOUBLE
+#define DOUBLE
 
 #ifdef DOUBLE
 	#define Real double
-	
-	#define ZERO 0.0
+
+ 	#define ZERO 0.0
 	#define ONE 1.0
 	#define TWO 2.0
-	#define FOUR 4.0
+ 	#define FOUR 4.0
+
+ 	#define SMALL 1.0e-10;
+
+	/** Reynolds number */
+	const Real Re_num = 1000.0;
+
+	/** SOR relaxation parameter */
+	const Real omega = 1.7;
+
+	/** Discretization mixture parameter (gamma) */
+	const Real mix_param = 0.9;
+
+	/** Safety factor for time step modification */
+	const Real tau = 0.5;
+
+	/** Body forces in x- and y- directions */
+	const Real gx = 0.0;
+	const Real gy = 0.0;
+
+	/** Domain size (non-dimensional) */
+	#define xLength 1.0
+	#define yLength 1.0
 #else
 	#define Real float
+
 	// replace double functions with float versions
+	#undef fmin
 	#define fmin fminf
+	#undef fmax
 	#define fmax fmaxf
+	#undef fabs
 	#define fabs fabsf
-	
-	#define ZERO 0.0f
+	#undef sqrt
+	#define sqrt sqrtf
+
+ 	#define ZERO 0.0f
 	#define ONE 1.0f
 	#define TWO 2.0f
 	#define FOUR 4.0f
+	#define SMALL 1.0e-10f;
+
+		/** Reynolds number */
+	const Real Re_num = 1000.0f;
+
+	/** SOR relaxation parameter */
+	const Real omega = 1.7f;
+
+	/** Discretization mixture parameter (gamma) */
+	const Real mix_param = 0.9f;
+
+	/** Safety factor for time step modification */
+	const Real tau = 0.5f;
+
+	/** Body forces in x- and y- directions */
+	const Real gx = 0.0f;
+	const Real gy = 0.0f;
+
+	/** Domain size (non-dimensional) */
+	#define xLength 1.0f
+	#define yLength 1.0f
 #endif
 
-/** Use shared memory */
-#define SHARED
-
-/** Use atomic operations to calculate residual, only for SINGLE PRECISION */
-//#define ATOMIC
-
-#if defined (ATOMIC) && defined (DOUBLE)
-# error double precision atomic operations not supported
-#endif
-
-typedef unsigned int uint;
-
-/** Reynolds number */
-#ifdef DOUBLE
-const Real Re_num = 1000.0;
-#else
-const Real Re_num = 1000.0f;
-#endif
-
-/** SOR relaxation parameter */
-#ifdef DOUBLE
-const Real omega = 1.7;
-#else
-const Real omega = 1.7f;
-#endif
-
-/** Discretization mixture parameter */
-#ifdef DOUBLE
-const Real mix_param = 0.9;
-#else
-const Real mix_param = 0.9f;
-#endif
-
-/** Safety factor for time step modification */
-#ifdef DOUBLE
-const Real tau = 0.5;
-#else
-const Real tau = 0.5f;
-#endif
-
-/** Body forces in x- and y- directions */
-const Real gx = ZERO;
-const Real gy = ZERO;
-
-/** Normalized lid velocity */
-const Real uStar = ONE;
-
-/** Domain size (non-dimensional) */
-#define xLength ONE
-#define yLength ONE
 
 /** Mesh sizes */
 const Real dx = xLength / NUM;
 const Real dy = yLength / NUM;
 
 /** Max macro (type safe, from GNU) */
-#define MAX(a,b) ({ __typeof__ (a) _a = (a); __typeof__ (b) _b = (b); _a > _b ? _a : _b; })
+//#define MAX(a,b) ({ __typeof__ (a) _a = (a); __typeof__ (b) _b = (b); _a > _b ? _a : _b; })
 
 /** Min macro (type safe) */
-#define MIN(a,b) ({ __typeof__ (a) _a = (a); __typeof__ (b) _b = (b); _a < _b ? _a : _b; })
-
-#ifdef ATOMIC
-// need bijection between float and unsigned int to use atomicMax()
-
-static __inline__ __device__ unsigned int floatFlip (float theFloat)
-{
-	unsigned int mask = (__float_as_int(theFloat) >> 31) | 0x80000000;
-	return __float_as_int(theFloat) ˆ mask;
-}
-
-inline __host__ float invFloatFlip (unsigned int theUint)
-{
-	unsigned int mask = ((theUint >> 31) - 1) | 0x80000000;
-	return __int_as_float((int)(theUint ^ mask));
-}
-#endif
+//#define MIN(a,b) ({ __typeof__ (a) _a = (a); __typeof__ (b) _b = (b); _a < _b ? _a : _b; })
+
+// map two-dimensional indices to one-dimensional memory
+#define u(I, J) u[((I) * ((NUM) + 2)) + (J)]
+#define v(I, J) v[((I) * ((NUM) + 2)) + (J)]
+#define F(I, J) F[((I) * ((NUM) + 2)) + (J)]
+#define G(I, J) G[((I) * ((NUM) + 2)) + (J)]
+#define pres_red(I, J) pres_red[((I) * ((NUM_2) + 2)) + (J)]
+#define pres_black(I, J) pres_black[((I) * ((NUM_2) + 2)) + (J)]
+
+///////////////////////////////////////////////////////////////////////////////
+__host__
+void set_BCs_host (Real* u, Real* v) 
+{
+	int ind;
+
+	// loop through rows and columns
+	for (ind = 0; ind < NUM + 2; ++ind) {
+
+		// left boundary
+		u(0, ind) = ZERO;
+		v(0, ind) = -v(1, ind);
+
+		// right boundary
+		u(NUM, ind) = ZERO;
+		v(NUM + 1, ind) = -v(NUM, ind);
+
+		// bottom boundary
+		u(ind, 0) = -u(ind, 1);
+		v(ind, 0) = ZERO;
+
+		// top boundary
+		u(ind, NUM + 1) = TWO - u(ind, NUM);
+		v(ind, NUM) = ZERO;
+
+		if (ind == NUM) {
+			// left boundary
+			u(0, 0) = ZERO;
+			v(0, 0) = -v(1, 0);
+			u(0, NUM + 1) = ZERO;
+			v(0, NUM + 1) = -v(1, NUM + 1);
+
+			// right boundary
+			u(NUM, 0) = ZERO;
+			v(NUM + 1, 0) = -v(NUM, 0);
+			u(NUM, NUM + 1) = ZERO;
+			v(NUM + 1, NUM + 1) = -v(NUM, NUM + 1);
+
+			// bottom boundary
+			u(0, 0) = -u(0, 1);
+			v(0, 0) = ZERO;
+			u(NUM + 1, 0) = -u(NUM + 1, 1);
+			v(NUM + 1, 0) = ZERO;
+
+			// top boundary
+			u(0, NUM + 1) = TWO - u(0, NUM);
+			v(0, NUM) = ZERO;
+			u(NUM + 1, NUM + 1) = TWO - u(NUM + 1, NUM);
+			v(ind, NUM + 1) = ZERO;
+		} // end if
+
+	} // end for
+
+} // end set_BCs_host
+
+///////////////////////////////////////////////////////////////////////////////
+__global__
+void set_BCs (Real* u, Real* v) 
+{
+	int ind = (blockIdx.x * blockDim.x) + threadIdx.x + 1;
+
+	// left boundary
+	u(0, ind) = ZERO;
+	v(0, ind) = -v(1, ind);
+
+	// right boundary
+	u(NUM, ind) = ZERO;
+	v(NUM + 1, ind) = -v(NUM, ind);
+
+	// bottom boundary
+	u(ind, 0) = -u(ind, 1);
+	v(ind, 0) = ZERO;
+
+	// top boundary
+	u(ind, NUM + 1) = TWO - u(ind, NUM);
+	v(ind, NUM) = ZERO;
+
+	if (ind == NUM) {
+		// left boundary
+		u(0, 0) = ZERO;
+		v(0, 0) = -v(1, 0);
+		u(0, NUM + 1) = ZERO;
+		v(0, NUM + 1) = -v(1, NUM + 1);
+
+		// right boundary
+		u(NUM, 0) = ZERO;
+		v(NUM + 1, 0) = -v(NUM, 0);
+		u(NUM, NUM + 1) = ZERO;
+		v(NUM + 1, NUM + 1) = -v(NUM, NUM + 1);
+
+		// bottom boundary
+		u(0, 0) = -u(0, 1);
+		v(0, 0) = ZERO;
+		u(NUM + 1, 0) = -u(NUM + 1, 1);
+		v(NUM + 1, 0) = ZERO;
+
+		// top boundary
+		u(0, NUM + 1) = TWO - u(0, NUM);
+		v(0, NUM) = ZERO;
+		u(NUM + 1, NUM + 1) = TWO - u(NUM + 1, NUM);
+		v(ind, NUM + 1) = ZERO;
+	} // end if
+
+} // end set_BCs
 
 ///////////////////////////////////////////////////////////////////////////////
 
 __global__ 
-void calculate_F (const Real * u, const Real * v, const Real dt, 
-									 Real * F)
+void calculate_F (const Real dt, const Real* u, const Real* v,
+				  Real* F) 
 {	
-	int row = (blockIdx.y * blockDim.y) + threadIdx.y;
-	int col = (blockIdx.x * blockDim.x) + threadIdx.x;
-	
-	if (col == (NUM - 1)) {
+	int row = (blockIdx.x * blockDim.x) + threadIdx.x + 1;
+	int col = (blockIdx.y * blockDim.y) + threadIdx.y + 1;
+	
+	if (col == NUM) {
 		// right boundary, F_ij = u_ij
-		F[((NUM - 1) * NUM) + row] = u[((NUM - 1) * NUM) + row];
+		// also do left boundary
+		F(0, row) = u(0, row);
+		F(NUM, row) = u(NUM, row);
 	} else {
 		
-		// u and v velocities
-		Real u_ij = u[(col * NUM) + row];
-		Real u_ip1j = u[((col + 1) * NUM) + row];
-		
-		Real v_ij = v[(col * NUM) + row];
-		Real v_ip1j = v[((col + 1) * NUM) + row];
-		
-		// left boundary
-		Real u_im1j;
-		if (col == 0) {
-			u_im1j = ZERO;
-		} else {
-			u_im1j = u[((col - 1) * NUM) + row];
-		}
-		
-		// bottom boundary
-		Real u_ijm1, v_ijm1, v_ip1jm1;
-		if (row == 0) {
-			u_ijm1 = -u_ij;
-			v_ijm1 = ZERO;
-			v_ip1jm1 = ZERO;
-		} else {
-			u_ijm1 = u[(col * NUM) + row - 1];
-			v_ijm1 = v[(col * NUM) + row - 1];
-			v_ip1jm1 = v[((col + 1) * NUM) + row - 1];
-		}
-		
-		// top boundary
-		Real u_ijp1;
-		if (row == (NUM - 1)) {
-			u_ijp1 = (TWO * uStar) - u_ij;
-		} else {
-			u_ijp1 = u[(col * NUM) + row + 1];
-		}
+		// u velocities
+		Real u_ij = u(col, row);
+		Real u_ip1j = u(col + 1, row);
+		Real u_ijp1 = u(col, row + 1);
+		Real u_im1j = u(col - 1, row);
+		Real u_ijm1 = u(col, row - 1);
+
+		// v velocities
+		Real v_ij = v(col, row);
+		Real v_ip1j = v(col + 1, row);
+		Real v_ijm1 = v(col, row - 1);
+		Real v_ip1jm1 = v(col + 1, row - 1);
 		
 		// finite differences
 		Real du2dx, duvdy, d2udx2, d2udy2;
 
 		du2dx = (((u_ij + u_ip1j) * (u_ij + u_ip1j) - (u_im1j + u_ij) * (u_im1j + u_ij))
-						+ mix_param * (fabs(u_ij + u_ip1j) * (u_ij - u_ip1j)
-						- fabs(u_im1j + u_ij) * (u_im1j - u_ij))) / (FOUR * dx);
+				+ mix_param * (fabs(u_ij + u_ip1j) * (u_ij - u_ip1j)
+				- fabs(u_im1j + u_ij) * (u_im1j - u_ij))) / (FOUR * dx);
 		duvdy = ((v_ij + v_ip1j) * (u_ij + u_ijp1) - (v_ijm1 + v_ip1jm1) * (u_ijm1 + u_ij)
-					+ mix_param * (fabs(v_ij + v_ip1j) * (u_ij - u_ijp1)
-					- fabs(v_ijm1 + v_ip1jm1) * (u_ijm1 - u_ij))) / (FOUR * dy);
-	 	d2udx2 = (u_ip1j - (TWO * u_ij) + u_im1j) / (dx * dx);
-	  	d2udy2 = (u_ijp1 - (TWO * u_ij) + u_ijm1) / (dy * dy);
-
-		F[(col * NUM) + row] = u_ij + dt * (((d2udx2 + d2udy2) / Re_num) - du2dx - duvdy + gx);
+				+ mix_param * (fabs(v_ij + v_ip1j) * (u_ij - u_ijp1)
+				- fabs(v_ijm1 + v_ip1jm1) * (u_ijm1 - u_ij))) / (FOUR * dy);
+		d2udx2 = (u_ip1j - (TWO * u_ij) + u_im1j) / (dx * dx);
+		d2udy2 = (u_ijp1 - (TWO * u_ij) + u_ijm1) / (dy * dy);
+
+		F(col, row) = u_ij + dt * (((d2udx2 + d2udy2) / Re_num) - du2dx - duvdy + gx);
 		
 	} // end if
 		

///////////////////////////////////////////////////////////////////////////////

__global__ 
void calculate_G (const Real * u, const Real * v, const Real dt, 
									 Real * G)
{
	int row = (blockIdx.y * blockDim.y) + threadIdx.y;
	int col = (blockIdx.x * blockDim.x) + threadIdx.x;
	
	if (row == (NUM - 1)) {
		G[(col * NUM) + NUM - 1] = v[(col * NUM) + NUM - 1];
	} else {
		
		// u and v velocities
		Real u_ij = u[(col * NUM) + row];
		Real u_ijp1 = u[(col * NUM) + row + 1];
		
		Real v_ij = v[(col * NUM) + row];
		Real v_ijp1 = v[(col * NUM) + row + 1];
		
		// bottom boundary
		Real v_ijm1;
		if (row == 0) {
			v_ijm1 = ZERO;
		} else {
			v_ijm1 = v[(col * NUM) + row - 1];
		}
		
		// left boundary
		Real v_im1j, u_im1j, u_im1jp1;
		if (col == 0) {
			v_im1j = -v_ij;
			u_im1j = ZERO;
			u_im1jp1 = ZERO;
		} else {
			v_im1j = v[((col - 1) * NUM) + row];
			u_im1j = u[((col - 1) * NUM) + row];
			u_im1jp1 = u[((col - 1) * NUM) + row + 1];
		}
		
		// right boundary
		Real v_ip1j;
		if (col == (NUM - 1)) {
			v_ip1j = -v_ij;
		} else {
			v_ip1j = v[((col + 1) * NUM) + row];
		}
		
		// finite differences
		Real dv2dy, duvdx, d2vdx2, d2vdy2;
	
		dv2dy = ((v_ij + v_ijp1) * (v_ij + v_ijp1) - (v_ijm1 + v_ij) * (v_ijm1 + v_ij)
		  		+ mix_param * (fabs(v_ij + v_ijp1) * (v_ij - v_ijp1)
					- fabs(v_ijm1 + v_ij) * (v_ijm1 - v_ij))) / (FOUR * dy);
		duvdx = ((u_ij + u_ijp1) * (v_ij + v_ip1j) - (u_im1j + u_im1jp1) * (v_im1j + v_ij)
					+ mix_param * (fabs(u_ij + u_ijp1) * (v_ij - v_ip1j) 
					- fabs(u_im1j + u_im1jp1) * (v_im1j - v_ij))) / (FOUR * dx);
	  	d2vdx2 = (v_ip1j - (TWO * v_ij) + v_im1j) / (dx * dx);
	  	d2vdy2 = (v_ijp1 - (TWO * v_ij) + v_ijm1) / (dy * dy);

		G[(col * NUM) + row] = v_ij + dt * (((d2vdx2 + d2vdy2) / Re_num) - dv2dy - duvdx + gy);
			
	} // end if
		

int main (void)
{
	// iterations for Red-Black Gauss-Seidel with SOR
	int iter = 0;
	const int it_max = 10000;
	
	// SOR iteration tolerance
	const Real tol = 0.001;
