Fixed HLL

SWESimulators/Common.py

@@ -53,9 +53,9 @@ class CUDAArray2D:
         self.ny_halo = ny + 2*halo_y
         
         #Make sure data is in proper format
-        assert(np.issubdtype(data.dtype, np.float32))
+        assert(np.issubdtype(data.dtype, np.float32), "Wrong datatype: %s" % str(data.dtype))
-        assert(not np.isfortran(data))
+        assert(not np.isfortran(data), "Wrong datatype (Fortran, expected C)")
-        assert(data.shape == (self.ny_halo, self.nx_halo))
+        assert(data.shape == (self.ny_halo, self.nx_halo), "Wrong data shape: %s" % str(data.shape))
 
         #Upload data to the device
         self.data = pycuda.gpuarray.to_gpu_async(data, stream=stream)

SWESimulators/FORCE.py

@@ -94,8 +94,8 @@ class FORCE:
         #Compute kernel launch parameters
         self.local_size = (block_width, block_height, 1) 
         self.global_size = ( \
-                       int(np.ceil(self.nx / float(self.local_size[0])) * self.local_size[0]), \
+                       int(np.ceil(self.nx / float(self.local_size[0]))), \
-                       int(np.ceil(self.ny / float(self.local_size[1])) * self.local_size[1]) \
+                       int(np.ceil(self.ny / float(self.local_size[1]))) \
                       ) 
     
     

SWESimulators/FORCE_kernel.cu

@@ -110,18 +110,6 @@ __global__ void FORCEKernel(
         float* hu1_ptr_, int hu1_pitch_,
         float* hv1_ptr_, int hv1_pitch_) {
     
-    //Index of thread within block
-    const int tx = get_local_id(0);
-    const int ty = get_local_id(1);
-    
-    //Index of block within domain
-    const int bx = get_local_size(0) * get_group_id(0);
-    const int by = get_local_size(1) * get_group_id(1);
-
-    //Index of cell within domain
-    const int ti = get_global_id(0) + 1; //Skip global ghost cells, i.e., +1
-    const int tj = get_global_id(1) + 1;
-    
     __shared__ float Q[3][block_height+2][block_width+2];
     __shared__ float F[3][block_height+1][block_width+1];
     

SWESimulators/HLL.py

@@ -34,7 +34,6 @@ from SWESimulators import Common
 
 
 
-
 """
 Class that solves the SW equations using the Harten-Lax -van Leer approximate Riemann solver
 """
@@ -43,8 +42,8 @@ class HLL:
     """
     Initialization routine
     h0: Water depth incl ghost cells, (nx+1)*(ny+1) cells
-    u0: Initial momentum along x-axis incl ghost cells, (nx+1)*(ny+1) cells
+    hu0: Initial momentum along x-axis incl ghost cells, (nx+1)*(ny+1) cells
-    v0: Initial momentum along y-axis incl ghost cells, (nx+1)*(ny+1) cells
+    hv0: Initial momentum along y-axis incl ghost cells, (nx+1)*(ny+1) cells
     nx: Number of cells along x-axis
     ny: Number of cells along y-axis
     dx: Grid cell spacing along x-axis (20 000 m)
@@ -90,8 +89,8 @@ class HLL:
         #Compute kernel launch parameters
         self.local_size = (block_width, block_height, 1) 
         self.global_size = ( \
-                       int(np.ceil(self.nx / float(self.local_size[0])) * self.local_size[0]), \
+                       int(np.ceil(self.nx / float(self.local_size[0]))), \
-                       int(np.ceil(self.ny / float(self.local_size[1])) * self.local_size[1]) \
+                       int(np.ceil(self.ny / float(self.local_size[1]))) \
                       ) 
     
     

SWESimulators/HLL_kernel.cu

@@ -59,7 +59,7 @@ void computeFluxF(float Q[3][block_height+2][block_width+2],
 
 
 /**
-  * Computes the flux along the x axis for all faces
+  * Computes the flux along the y axis for all faces
   */
 __device__
 void computeFluxG(float Q[3][block_height+2][block_width+2],
@@ -148,6 +148,8 @@ __global__ void HLLKernel(
     __syncthreads();
     
     
+    //Q[0][get_local_id(1) + 1][get_local_id(0) + 1] += 0.1;
+    
     
     
     // Write to main memory for all internal cells

SWESimulators/LxF.py

@@ -35,12 +35,8 @@ from SWESimulators import Common
 
 
 
-        
-        
-
-
 """
-Class that solves the SW equations using the Forward-Backward linear scheme
+Class that solves the SW equations using the Lax Friedrichs scheme
 """
 class LxF:
 
@@ -63,7 +59,7 @@ class LxF:
                  g, \
                  block_width=16, block_height=16):
         #Create a CUDA stream
-        self.stream = cuda.Stream()
+        self.stream = None #cuda.Stream()
 
         #Get kernels
         self.lxf_module = Common.get_kernel("LxF_kernel.cu", block_width, block_height)
@@ -94,8 +90,8 @@ class LxF:
         #Compute kernel launch parameters
         self.local_size = (block_width, block_height, 1) 
         self.global_size = ( \
-                       int(np.ceil(self.nx / float(self.local_size[0])) * self.local_size[0]), \
+                       int(np.ceil(self.nx / float(self.local_size[0]))), \
-                       int(np.ceil(self.ny / float(self.local_size[1])) * self.local_size[1]) \
+                       int(np.ceil(self.ny / float(self.local_size[1]))) \
                       ) 
     
     

SWESimulators/common.cu

@@ -90,7 +90,6 @@ inline __device__ float3 operator+(const float3 a, const float3 b) {
     return make_float3(a.x+b.x, a.y+b.y, a.z+b.z);
 }
 
-
 inline __device__ __host__ float clamp(const float f, const float a, const float b) {
     return fmaxf(a, fminf(f, b));
 }
@@ -834,13 +833,13 @@ __device__ float3 WAF_1D_flux(const float3 Q_l2, const float3 Q_l1, const float3
     
     // Compute the r parameters for the flux limiter
     const float rh_1 = (c_1 > 0.0f) ? rh_m : rh_p; 
-    const float rv_1 = (c_1 > 0.0f) ? rv_m : rv_p; 
+    //const float rv_1 = (c_1 > 0.0f) ? rv_m : rv_p; 
     
-    const float rh_2 = (c_2 > 0.0f) ? rh_m : rh_p; 
+    //const float rh_2 = (c_2 > 0.0f) ? rh_m : rh_p; 
     const float rv_2 = (c_2 > 0.0f) ? rv_m : rv_p; 
     
     const float rh_3 = (c_3 > 0.0f) ? rh_m : rh_p;
-    const float rv_3 = (c_3 > 0.0f) ? rv_m : rv_p;
+    //const float rv_3 = (c_3 > 0.0f) ? rv_m : rv_p;
     
     // Compute the limiter
     // We use h for the nonlinear waves, and v for the middle shear wave