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


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

/** Function to update pressure for red cells
 * 
 * \param[in]			dt					time-step size
 * \param[in]			F						array of discretized x-momentum eqn terms
 * \param[in]			G						array of discretized y-momentum eqn terms
 * \param[in]			pres_black	pressure values of black cells
 * \param[inout]	pres_red		pressure values of red cells
 * \param[inout]	norm_L2			variable holding summed residuals
 */
__global__ 
void red_kernel (const Real dt, const Real * F, const Real * G, const Real * pres_black,
								 Real * pres_red, Real * norm_L2)
{
	int row = (blockIdx.y * blockDim.y) + threadIdx.y;
	int col = (blockIdx.x * blockDim.x) + threadIdx.x;
	
	int ind_red = (col * (NUM >> 1)) + row;  					// local (red) index
	int ind = (col * NUM) + (2 * row) + (col & 1);		// global index

	Real p_ij = pres_red[ind_red];			

	// left boundary
	Real p_im1j;
	Real F_im1j;
	if (col == 0) {
		p_im1j = p_ij;
		F_im1j = ZERO;
	} else {
		p_im1j = pres_black[((col - 1) * (NUM >> 1)) + row];
		F_im1j = F[((col - 1) * NUM) + (2 * row) + (col & 1)];
	}

	// right boundary
	Real p_ip1j;
	if (col == (NUM - 1)) {
		p_ip1j = p_ij;
	} else {
		p_ip1j = pres_black[((col + 1) * (NUM >> 1)) + row];
	}

	// bottom boundary
	Real p_ijm1;
	Real G_ijm1;
	if (((2 * row) + (col & 1)) == 0) {
		p_ijm1 = p_ij;
		G_ijm1 = ZERO;
	} else {
		p_ijm1 = pres_black[(col * (NUM >> 1)) + row - ((col + 1) & 1)];
		G_ijm1 = G[(col * NUM) + (2 * row) + (col & 1) - 1];
	}

	// top boundary
	Real p_ijp1;
	if (((2 * row) + (col & 1)) == (NUM - 1)) {
		p_ijp1 = p_ij;
	} else {
		p_ijp1 = pres_black[(col * (NUM >> 1)) + row + (col & 1)];
	}

	// right-hand side
	Real rhs = (((F[ind] - F_im1j) / dx) + ((G[ind] - G_ijm1) / dy)) / dt;

	pres_red[ind_red] = p_ij * (ONE - omega) + omega * (
										  ((p_ip1j + p_im1j) / (dx * dx)) + ((p_ijp1 + p_ijm1) / (dy * dy)) - 
										  rhs) / ((TWO / (dx * dx)) + (TWO / (dy * dy)));
	
	// calculate residual (reuse rhs variable)
	rhs = ((p_ip1j - (TWO * p_ij) + p_im1j) / (dx * dx)) + ((p_ijp1 - (TWO * p_ij) + p_ijm1) / (dy * dy)) - rhs;
	
	#ifdef SHARED
		// store residual for block
		__shared__ Real res_cache[BLOCK_SIZE];
	
		// store squared residual from each thread
		res_cache[threadIdx.y] = rhs * rhs;
	
		// synchronize threads in block
		__syncthreads();
	
		// add up squared residuals for block
		int i = BLOCK_SIZE >> 1;
		while (i != 0) {
			if (threadIdx.y < i) {
				res_cache[threadIdx.y] += res_cache[threadIdx.y + i];
			}
			__syncthreads();
			i >>= 1;
		}
	
		// store block's summed residuals
		if (threadIdx.y == 0) {
			#ifdef ATOMIC
				atomicAdd (norm_L2, res_cache[0]);
			#else
				norm_L2[blockIdx.x + (gridDim.x * blockIdx.y)] = res_cache[0];
			#endif
		}
	#else
		norm_L2[ind_red] = rhs * rhs;
	#endif
	
} // end red_kernel

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

/** Function to update pressure for black cells
 * 
 * \param[in]			dt					time-step size
 * \param[in]			F						array of discretized x-momentum eqn terms
 * \param[in]			G						array of discretized y-momentum eqn terms
 * \param[in]			pres_red		pressure values of red cells
 * \param[inout]	pres_black	pressure values of black cells
 * \param[inout]	norm_L2			variable holding summed residuals
 */
__global__ 
void black_kernel (const Real dt, const Real * F, const Real * G, const Real * pres_red,
									 Real * pres_black, Real * norm_L2)
{
	int row = (blockIdx.y * blockDim.y) + threadIdx.y;
	int col = (blockIdx.x * blockDim.x) + threadIdx.x;
	
	int ind_black = (col * (NUM >> 1)) + row;  						// local (black) index
	int ind = (col * NUM) + (2 * row) + ((col + 1) & 1);	// global index
	
	Real p_ij = pres_black[ind_black];
	
	// left boundary
	Real p_im1j;
	Real F_im1j;
	if (col == 0) {
		p_im1j = p_ij;
		F_im1j = ZERO;
	} else {
		p_im1j = pres_red[((col - 1) * (NUM >> 1)) + row];
		F_im1j = F[((col - 1) * NUM) + (2 * row) + ((col + 1) & 1)];
	}
	
	// right boundary
	Real p_ip1j;
	if (col == (NUM - 1)) {
		p_ip1j = p_ij;
	} else {
		p_ip1j = pres_red[((col + 1) * (NUM >> 1)) + row];
	}
	
	// bottom boundary
	Real p_ijm1;
	Real G_ijm1;
	if (((2 * row) + ((col + 1) & 1)) == 0) {
		p_ijm1 = p_ij;
		G_ijm1 = ZERO;
	} else {
		p_ijm1 = pres_red[(col * (NUM >> 1)) + row - (col & 1)];
		G_ijm1 = G[(col * NUM) + (2 * row) + ((col + 1) & 1) - 1];
	}
	
	// top boundary
	Real p_ijp1;
	if (((2 * row) + ((col + 1) & 1)) == (NUM - 1)) {
		p_ijp1 = p_ij;
	} else {
		p_ijp1 = pres_red[(col * (NUM >> 1)) + row + ((col + 1) & 1)];
	}
	
	// right-hand side
	Real rhs = (((F[ind] - F_im1j) / dx) + ((G[ind] - G_ijm1) / dy)) / dt;
	
	pres_black[ind_black] = p_ij * (ONE - omega) + omega * (
										  		((p_ip1j + p_im1j) / (dx * dx)) + ((p_ijp1 + p_ijm1) / (dy * dy)) - 
										  		rhs) / ((TWO / (dx * dx)) + (TWO / (dy * dy)));
	
	// calculate residual (reuse rhs variable)
	rhs = ((p_ip1j - (TWO * p_ij) + p_im1j) / (dx * dx)) + ((p_ijp1 - (TWO * p_ij) + p_ijm1) / (dy * dy)) - rhs;
	
	#ifdef SHARED
		// store residual for block
		__shared__ Real res_cache[BLOCK_SIZE];
	
		// store squared residual from each thread
		res_cache[threadIdx.y] = rhs * rhs;
	
		// synchronize threads in block
		__syncthreads();
	
		// add up squared residuals for block
		int i = BLOCK_SIZE >> 1;
		while (i != 0) {
			if (threadIdx.y < i) {
				res_cache[threadIdx.y] += res_cache[threadIdx.y + i];
			}
			__syncthreads();
			i >>= 1;
		}
	
		// store block's summed residuals
		if (threadIdx.y == 0) {
			#ifdef ATOMIC
				atomicAdd (norm_L2, res_cache[0]);
			#else
				norm_L2[blockIdx.x + (gridDim.x * blockIdx.y)] = res_cache[0];
			#endif
		}
	#else
		norm_L2[ind_black] = rhs * rhs;
	#endif
	
} // end black_kernel

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

__global__ 
void calculate_u (const Real dt, const Real * F, 
									const Real * pres_red, const Real * pres_black, 
									#ifdef ATOMIC
									Real * u, unsigned int * max_u_d)
									#else
									Real * u, Real * max_u_d)
									#endif
{
	int row = (blockIdx.y * blockDim.y) + threadIdx.y;
	int col = (blockIdx.x * blockDim.x) + threadIdx.x;
	
	#ifdef SHARED
		#ifdef ATOMIC
			__shared__ unsigned int res_cache[BLOCK_SIZE];
			res_cache[threadIdx.y] = 0;
		#else
			__shared__ Real res_cache[BLOCK_SIZE];
			res_cache[threadIdx.y] = ZERO;
		#endif
	#endif
	
	if (col != (NUM - 1)) {
		int ind = (col * NUM) + row;
		
		Real p_ij, p_ip1j;
		if (((row + col) & 1) == 0) {
			// red pressure cell
			p_ij = pres_red[(col * (NUM >> 1)) + ((row - (col & 1)) >> 1)];
			
			// p_ip1j is black cell
			p_ip1j = pres_black[((col + 1) * (NUM >> 1)) + ((row - (col & 1)) >> 1)];
		} else {
			// black pressure cell
			p_ij = pres_black[(col * (NUM >> 1)) + ((row - ((col + 1) & 1)) >> 1)];
			
			// p_ip1j is red cell
			p_ip1j = pres_red[((col + 1) * (NUM >> 1)) + ((row - ((col + 1) & 1)) >> 1)];
		}
		
		//u[ind] = F[ind] - (dt * (p_ip1j - p_ij) / dx);
		Real u_ij = F[ind] - (dt * (p_ip1j - p_ij) / dx);
		
		u[ind] = u_ij;
		
		#ifdef SHARED
		// store maximum u for block from each thread
			#ifdef ATOMIC
				res_cache[threadIdx.y] = floatFlip (fabs(u_ij));
			#else
				res_cache[threadIdx.y] = fabs(u_ij);
			#endif
		
		// synchronize threads in block
		__syncthreads();

		// add up squared residuals for block
		int i = BLOCK_SIZE >> 1;
		while (i != 0) {
			if (threadIdx.y < i) {
				res_cache[threadIdx.y] = fmax(res_cache[threadIdx.y], res_cache[threadIdx.y + i]);
			}
			__syncthreads();
			i >>= 1;
		}

		// store block's summed residuals
		if (threadIdx.y == 0) {
			#ifdef ATOMIC
				atomicMax (max_u_d, res_cache[0]);
			#else
				max_u_d[blockIdx.x + (gridDim.x * blockIdx.y)] = res_cache[0];
			#endif
		}
		#endif
	} // end if
	
} // end calculate_u

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

__global__ 
void calculate_v (const Real dt, const Real * G, 
									const Real * pres_red, const Real * pres_black, 
									#ifdef ATOMIC
									Real * v, unsigned int * max_v_d)
									#else
									Real * v, Real * max_v_d)
									#endif
{
	int row = (blockIdx.y * blockDim.y) + threadIdx.y;
	int col = (blockIdx.x * blockDim.x) + threadIdx.x;
	
	#ifdef SHARED
		#ifdef ATOMIC
			__shared__ unsigned int res_cache[BLOCK_SIZE];
			res_cache[threadIdx.y] = 0;
		#else
			__shared__ Real res_cache[BLOCK_SIZE];
			res_cache[threadIdx.y] = ZERO;
		#endif
	#endif
	
	if (row != (NUM - 1)) {
		int ind = (col * NUM) + row;
		
		Real p_ij, p_ijp1;
		if (((row + col) & 1) == 0) {
			// red pressure cell
			p_ij = pres_red[(col * (NUM >> 1)) + ((row - (col & 1)) >> 1)];
			
			// p_ijp1 is black cell
			p_ijp1 = pres_black[(col * (NUM >> 1)) + ((row + 1 - ((col + 1) & 1)) >> 1)];
		} else {
			// black pressure cell
			p_ij = pres_black[(col * (NUM >> 1)) + ((row - ((col + 1) & 1)) >> 1)];
			
			// p_ijp1 is red cell
			p_ijp1 = pres_red[(col * (NUM >> 1)) + ((row + 1 - (col & 1)) >> 1)];
		}
		
		//v[ind] = G[ind] - (dt * (p_ijp1 - p_ij) / dy);
		Real v_ij = G[ind] - (dt * (p_ijp1 - p_ij) / dy);
		
		v[ind] = v_ij;
		
		#ifdef SHARED
		// store maximum v for block for each thread
		
			#ifdef ATOMIC
				res_cache[threadIdx.y] = floatFlip (fabs(v_ij));
			#else
				res_cache[threadIdx.y] = fabs(v_ij);
			#endif
		
		// synchronize threads in block
		__syncthreads();

		// add up squared residuals for block
		int i = BLOCK_SIZE >> 1;
		while (i != 0) {
			if (threadIdx.y < i) {
				res_cache[threadIdx.y] = fmax(res_cache[threadIdx.y], res_cache[threadIdx.y + i]);
			}
			__syncthreads();
			i >>= 1;
		}

		// store block's summed residuals
		if (threadIdx.y == 0) {
			#ifdef ATOMIC
				atomicMax (max_v_d, res_cache[0]);
			#else
				max_v_d[blockIdx.x + (gridDim.x * blockIdx.y)] = res_cache[0];
			#endif
		}
		#endif
	} // end if
	
} // end calculate_v


	free(F);
	free(G);
	
	free(max_u_arr);
	free(max_v_arr);
	free(res);
	free(pres_sum);
	
	checkCudaErrors (cudaDeviceReset());
	
	return 0;
}
