From 8f90f21afaad6d69295bcb11475a12826c30c2ab Mon Sep 17 00:00:00 2001 From: ph Date: Tue, 28 Nov 2017 21:55:35 +0100 Subject: [PATCH] transport lbm model works --- lbm_cl.py | 28 ++++++++++++++++++---------- lbm_kernels.cl | 12 +++++++----- 2 files changed, 25 insertions(+), 15 deletions(-) diff --git a/lbm_cl.py b/lbm_cl.py index 1bf94bc..8b65d1d 100755 --- a/lbm_cl.py +++ b/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() diff --git a/lbm_kernels.cl b/lbm_kernels.cl index c2ccad2..bdfe88e 100644 --- a/lbm_kernels.cl +++ b/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]; } -- GitLab