transport lbm model works

lbm_cl.py

@@ -25,8 +25,8 @@ _n = 4 * _m
 
 
 # grid size
-_nx = 64
-_ny = 64
+_nx = 1024
+_ny = 1024
 
 Lx = 1.
 Ly = 1.
@@ -38,11 +38,12 @@ _dy = Ly / _ny
 vel = np.array([1., 1.])
 
 # lattice speed
-_vmax = 2.
+_vmax = 3.
 
 
 # time stepping
-_Tmax = 0.5 / _vmax
+_Tmax = 10. / _vmax
+#_Tmax = 0.
 cfl = 1
 
 _dt = cfl * _dx / _vmax
@@ -103,7 +104,7 @@ def solve_ocl(m = _m, n = _n, nx = _nx, ny = _ny,
     event.wait()
 
     # number of animation frames
-    nbplots = 10
+    nbplots = 100
     itermax = int(np.floor(Tmax / dt))
     iterplot = int(itermax / nbplots)
 
@@ -127,9 +128,10 @@ def solve_ocl(m = _m, n = _n, nx = _nx, ny = _ny,
         print("iter=",iter, " t=",t, "elapsed (s)=",elapsed)
         if iter % iterplot == 0 and animate:
             cl.enqueue_copy(queue, fn_cpu, fn_gpu).wait()
-            wplot = np.reshape(fn_cpu, (nx, ny))
+            wplot = np.reshape(fn_cpu, (n, nx, ny))
             plt.clf()
-            plt.imshow(wplot,vmin=0, vmax=1)
+            #plt.imshow(np.sum(wplot, axis = 0),vmin=0, vmax=1)
+            plt.imshow(np.sum(wplot, axis = 0))
             plt.gca().invert_yaxis()
             plt.colorbar()
             plt.pause(0.01)
@@ -137,17 +139,23 @@ def solve_ocl(m = _m, n = _n, nx = _nx, ny = _ny,
     # copy OpenCL data to CPU and return the results
     cl.enqueue_copy(queue, fn_cpu, fn_gpu).wait()
 
-    wplot_gpu = np.reshape(fn_cpu,(nx, ny))
+    wplot_gpu = np.reshape(fn_cpu,(n, nx, ny))
     return wplot_gpu
 
 # gpu solve
-wplot_gpu = solve_ocl()
+wplot_gpu = solve_ocl(animate = True)
 plt.clf()
-plt.imshow(wplot_gpu, vmin=0, vmax=1)
+plt.imshow(np.sum(wplot_gpu,axis=0), vmin=0, vmax=1)
 plt.gca().invert_yaxis()
 plt.colorbar()
 plt.show()
 
+# for iv in range(4):
+#     plt.imshow(wplot_gpu[iv,:,:])
+#     plt.gca().invert_yaxis()
+#     plt.colorbar()
+#     plt.show()
+
 
 # check difference
 # plt.clf()

lbm_kernels.cl

@@ -35,7 +35,7 @@ void w2f(const float *w, float *f){
     float vnorm[2] = {dir[d][0], dir[d][1]};
     flux_phy(w, vnorm, flux);   
     for(int iv = 0; iv < _M; iv++){
-      f[d * _M + iv] = w[iv] + flux[iv] / _LAMBDA;
+      f[d * _M + iv] = w[iv] / 4 + flux[iv] / 2 / _LAMBDA;
     }
   }
 }
@@ -46,7 +46,7 @@ void f2w(const float *f, float *w){
   for(int iv = 0; iv < _M; iv++) w[iv] = 0;
   for(int d = 0; d < 4; d++){
     for(int iv = 0; iv < _M; iv++){
-      w[iv] = f[d * _M + iv];
+      w[iv] += f[d * _M + iv];
     }
   }
 }
@@ -86,6 +86,7 @@ __kernel void init_sol(__global  float *fn){
   for(int ik = 0; ik < _N; ik++){
     int imem = i + j * _NX + ik * ngrid;
     fn[imem] =  fnow[ik];
+    //fn[imem] = j;
   }
 
 }
@@ -107,8 +108,8 @@ __kernel void time_step(__global  float *fn, __global float *fnp1){
   
   // periodic shift in each direction
   for(int d = 0; d < 4; d++){
-    int iR = (i - dir[d][0]) % _NX;
-    int jR = (j - dir[d][1]) % _NY;
+    int iR = (i - dir[d][0] + _NX) % _NX;
+    int jR = (j - dir[d][1] + _NY) % _NY;
     
     for(int iv = 0; iv < _M; iv++){
       int ik = d * _M + iv;
@@ -129,9 +130,10 @@ __kernel void time_step(__global  float *fn, __global float *fnp1){
     fnext[ik] = 2 * fnext[ik] - fnow[ik];
 
   // save
-  for(int ik = 0; ik < _M; ik++){
+  for(int ik = 0; ik < _N; ik++){
       int imem = i + j * _NX + ik * ngrid;
       fnp1[imem] = fnext[ik];
+      //fnp1[imem] = fnow[ik];
   }