up

mhd/vortex_kernels.cl

@@ -310,6 +310,7 @@ __kernel void time_step(__global const real *fn, __global real *fnp1) {
       // if (iv == ePsi) {
       // fnext[ik] = 1.0 * fnext[ik];
       fnext[ik] = om * fnext[ik] - (om - 1) * fnow[ik];
+      //fnext[ik] = fnow[ik];
       //} else {
       // fnext[ik] = 1.9 * fnext[ik] - 0.9 * fnow[ik];
       // fnext[ik] = 2.0 * fnext[ik] - 1.0 * fnow[ik];
@@ -326,6 +327,78 @@ __kernel void time_step(__global const real *fn, __global real *fnp1) {
   }
 }
 
+
+// one time step of the LBM scheme
+__kernel void time_shift(__global const real *fn, __global real *fnp1) {
+
+  int id = get_global_id(0);
+
+  int i = id % _NX;
+  int j = id / _NX;
+
+  int ngrid = _NX * _NY;
+
+  real fnow[_N];
+
+  // shift of values in domain
+  for (int d = 0; d < 4; d++) {
+    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;
+      int imem = iR + jR * _NX + ik * ngrid;
+      // #ifdef dirichlet_updown
+      //       // dirichlet condition on up and down borders
+      //       // (values of border cells are unchanged)
+      //       if ((j == 0) || (j == _NY - 1)) {
+      //         imem = i + j * _NX + ik * ngrid;
+      //       }
+      // #elif defined dirichlet_leftright
+      //       // dirichlet condition on left and right borders
+      //       // (values of border cells are unchanged)
+      //       if ((i == 0) || (i == _NX - 1)) {
+      //         imem = i + j * _NX + ik * ngrid;
+      //       }
+      // #elif defined dirichlet
+      //       // dirichlet condition on all borders
+      //       // (values of border cells are unchanged)
+      //       if ((i == 0) || (i == _NX - 1) || (j == 0) || (j == _NY - 1)) {
+      //         imem = i + j * _NX + ik * ngrid;
+      //       }
+      // #endif
+      fnow[ik] = fn[imem];
+    }
+  }
+
+  real fnext[_N];
+
+  // transfer
+  for (int iv = 0; iv < _M; iv++) {
+    for (int d = 0; d < 4; d++) {
+      int ik = d * _M + iv;
+      // if (iv == ePsi) {
+      // fnext[ik] = 1.0 * fnext[ik];
+      fnext[ik] = fnow[ik];
+      //fnext[ik] = fnow[ik];
+      //} else {
+      // fnext[ik] = 1.9 * fnext[ik] - 0.9 * fnow[ik];
+      // fnext[ik] = 2.0 * fnext[ik] - 1.0 * fnow[ik];
+      // fnext[ik] = 1.0 * fnext[ik] - 0.0 * fnow[ik];
+      //}
+    }
+  }
+
+  // save
+  for (int ik = 0; ik < _N; ik++) {
+    int imem = i + j * _NX + ik * ngrid;
+    fnp1[imem] = fnext[ik];
+    //fnp1[imem] = fnow[ik];
+  }
+}
+
+
+
 // numerical divergence of B
 __kernel void div_b(__global real *fn, __global real *div_b) {
 

mhd/vortex_lbm.py

@@ -199,7 +199,7 @@ def solve_ocl(m=_m, n=_n, nx=_nx, ny=_ny, Lx=_Lx, Ly=_Ly, Tmax=_Tmax,
 
         t += dt
         #event = prg.time_step(queue, (nx * ny, ), None, wn_gpu, wnp1_gpu)
-        event = prg.time_step(queue, (nx * ny, ), None, fn_gpu, fnp1_gpu)
+        event = prg.time_shift(queue, (nx * ny, ), None, fn_gpu, fnp1_gpu)
         #event = prg.time_step(queue, (nx * ny, ), None, wn_gpu, wnp1_gpu, wait_for = [event])
         event.wait()
         queue.finish()