first commit: D2Q9 code to solve lod driven cavity test case

without using any for-loops.

Author: maruthinh

.gitignore

@@ -0,0 +1,6 @@
+#ignoring some folder related to python probably gets
+generated by spyder IDE
+D2Q9/__pycache__/
+.idea
+__pycache__
+*.png

README.md

@@ -0,0 +1,10 @@
+D2Q9 LBM code to solve lid driven cavity problem. This code has zero for loops
+
+To run the code 
+
+```
+python3.8 lid_driven_cavity.py
+```
+Velocity contour plot obtained from the solver. Simulation is carried out 
+with 64 x 64 grid points. Solution has converged to 1e-10
+

bc.py

@@ -0,0 +1,87 @@
+from d2q9 import D2Q9
+
+
+class BoundaryConditions(D2Q9):
+
+    def __init__(self, f):
+        self.f = f
+        self.f[:, 0, :] = self.f[:, 1, :]
+        self.f[:, -1, :] = self.f[:, -2, :]
+        self.f[:, :, 0] = self.f[:, :, 1]
+        self.f[:, :, -1] = self.f[:, :, -2]
+        D2Q9.__init__(self)
+        
+    def periodic_in_x(self):
+        """
+        To apply periodic boundary condition in x
+        """
+        self.f[1, 0, 1:-1] = self.f[1, -1, 1:-1]
+        self.f[8, 0, 1:-1] = self.f[8, -1, 1:-1]
+        self.f[5, 0, 1:-1] = self.f[5, -1, 1:-1]
+
+        self.f[6, -1, 1:-1] = self.f[6, 0, 1:-1]
+        self.f[3, -1, 1:-1] = self.f[3, 0, 1:-1]
+        self.f[7, -1, 1:-1] = self.f[7, 0, 1:-1]
+
+    def bounce_back_right(self):
+        """
+        Apply bounce back on right wall
+        """
+        self.f[3, -2, 2:-2] = self.f[1, -1, 2:-2]
+        self.f[6, -2, 2:-2] = self.f[8, -1, 2:-2]
+        self.f[7, -2, 2:-2] = self.f[5, -1, 2:-2]
+
+    def bounce_back_left(self):
+        """
+        Apply bounce back on right wall
+        """
+        self.f[8, 1, 2:-2] = self.f[6, 0, 2:-2]
+        self.f[1, 1, 2:-2] = self.f[3, 0, 2:-2]
+        self.f[5, 1, 2:-2] = self.f[7, 0, 2:-2]
+
+    def bounce_back_bottom(self):
+        """
+        Apply bounce back on right wall
+        """
+        self.f[2, 2:-2, 1] = self.f[4, 2:-2, 0]
+        self.f[6, 2:-2, 1] = self.f[8, 2:-2, 0]
+        self.f[5, 2:-2, 1] = self.f[7, 2:-2, 0]
+
+    def moving_wall_bc_top(self, u_wall):
+        """
+        Apply moving wall bounce back on top wall
+        """
+        density=self.f.sum(axis=0)
+        self.f[4, 1:-1, -2] = self.f[2, 1:-1, -1]
+        self.f[7, 1:-1, -2] = self.f[5, 1:-1, -1] + 6.0 * density[1:-1, -1] * self.w[7] * self.c[0][7] * u_wall
+        self.f[8, 1:-1, -2] = self.f[6, 1:-1, -1] + 6.0 * density[1:-1, -1] * self.w[8] * self.c[0][8] * u_wall
+
+    def bot_left_corner_correction(self):
+        """
+        Bottom left corner
+        """
+        self.f[2, 1, 1] = self.f[4, 1, 0]
+        self.f[5, 1, 1] = self.f[7, 1, 0]
+        self.f[1, 1, 1] = self.f[3, 1, 0]
+        self.f[8, 1, 1] = self.f[6, 1, 0]
+        self.f[6, 1, 1] = self.f[8, 1, 0]
+
+    def bot_right_corner_correction(self):
+        """
+        Bottom right corner
+        """
+        self.f[2, -2, 1] = self.f[4, -2, 0]
+        self.f[6, -2, 1] = self.f[8, -2, 0]
+        self.f[3, -2, 1] = self.f[1, -2, 0]
+        self.f[5, -2, 1] = self.f[7, -2, 0]
+        self.f[7, -2, 1] = self.f[5, -2, 0]
+
+    def top_left_corner_correction(self):
+        """
+        Top left corner
+        """
+        self.f[5, 1, -2] = self.f[7, 1, -1]
+        self.f[7, 1, -2] = self.f[5, 1, -1]
+        self.f[1, 1, -2] = self.f[3, 1, -1]
+        self.f[8, 1, -2] = self.f[6, 1, -1]
+        self.f[4, 1, -2] = self.f[2, 1, -1]

d2q9.py

@@ -0,0 +1,65 @@
+import numpy as np
+
+
+class D2Q9:
+
+    def __init__(self):
+        self.c = np.array([[0.0, 1.0, 0.0, -1.0, 0.0, 1.0, -1.0, -1.0, 1.0],
+                           [0.0, 0.0, 1.0, 0.0, -1.0, 1.0, 1.0, -1.0, -1.0]])
+        self.w = np.array([4.0 / 9.0, 1.0 / 9.0, 1.0 / 9.0, 1.0 / 9.0, 1.0 / 9.0,
+                           1.0 / 36.0, 1.0 / 36.0, 1.0 / 36.0, 1.0 / 36.0])
+        self.q = 9
+        self.T0 = 1 / 3.0
+        self.a = np.sqrt(self.T0)
+
+    def feq(self, moments):
+        """
+        To compute equilibrium distribution function.
+        """
+        
+        feq = np.zeros((self.q, moments.shape[1], moments.shape[2]))
+            
+        udotc = (np.array([[self.c[0]]]).T*moments[1]) + (np.array([[self.c[1]]]).T*moments[2])
+        usq = moments[1] * moments[1] + moments[2] * moments[2]
+        feq = (np.array([[self.w]]).T*moments[0]) * (1.0 + 3.0 * udotc - 1.5 * usq + 4.5 * udotc ** 2
+                                             + 4.5 * udotc ** 3 - 4.5 * usq * udotc)
+        return feq
+
+    def moments(self, f):
+        """
+        Compute moments (rho, u, v) from populations 
+        """
+        mom = np.zeros((3, f.shape[1], f.shape[2]))
+
+        mom[0] = f.sum(axis=0)
+        mom[1] = (np.array([[self.c[0]]]).T*f).sum(axis=0)/mom[0]
+        mom[2] = (np.array([[self.c[1]]]).T*f).sum(axis=0)/mom[0]
+
+        return mom
+
+    def collision(self, f, beta, alpha=2.0):
+        """
+        computes collision term 
+        """
+
+        feq = self.feq(self.moments(f))
+        
+        f += alpha * beta * (feq - f)
+        f.cumsum(axis=0)
+        
+
+    @staticmethod
+    def advection(f):
+        """
+        Streaming of velocities
+        """
+
+        f[3, :-1, :] = f[3, 1:, :]
+        f[4, :, :-1] = f[4, :, 1:]
+        f[7, :-1, :-1] = f[7, 1:, 1:]
+        f[8, 1:, :-1] = f[8, :-1, 1:]
+
+        f[1, -1:0:-1, ::-1] = f[1, -2::-1, ::-1]
+        f[2, ::-1, -1:0:-1] = f[2, ::-1, -2::-1]
+        f[5, -1:0:-1, -1:0:-1] = f[5, -2::-1, -2::-1]
+        f[6, -2:0:-1, -1:0:-1] = f[6, -1:1:-1, -2::-1]

lid_driven_cavity.py

@@ -0,0 +1,81 @@
+import numpy as np
+import matplotlib.pyplot as plt
+import time
+
+from d2q9 import D2Q9
+from bc import BoundaryConditions
+from sim_params import SimParams
+
+nx = 64
+ny = 64
+ghost_x = 2
+ghost_y = 2
+tot_nx = nx + ghost_x
+tot_ny = ny + ghost_y
+ndim = 2
+model = D2Q9()
+sim_params = SimParams(mach=0.1, reynolds_num=100, ref_len=1.0, res=0.02)
+q = model.q
+
+rho_ini = 1.0
+p = rho_ini * model.T0
+u_ini = 0.0
+v_ini = 0.0
+
+#initialization of arrays required
+f = np.zeros((q, tot_nx, tot_ny))
+moments = np.zeros((3, tot_nx, tot_ny))
+moments[0] = rho_ini
+moments[1] = u_ini
+moments[2] = v_ini
+
+u0 = np.zeros((tot_nx, tot_ny))
+v0 = np.zeros((tot_nx, tot_ny))
+
+
+f = model.feq(moments)
+
+iter = 0
+err = 1.0
+start = time.time()
+
+#main iteration loop
+while(err>1e-10):
+    iter+=1
+    model.collision(f, sim_params.beta)
+    bc = BoundaryConditions(f)
+    model.advection(f)
+    bc.bounce_back_left()
+    bc.bounce_back_bottom()
+    bc.bounce_back_right()
+    bc.moving_wall_bc_top(sim_params.ref_vel)
+
+    if iter % 1000 == 0:
+        moments=model.moments(f)
+        e1 = np.sum(np.sqrt((moments[1] - u0) ** 2 + (moments[2] - v0) ** 2))
+        e2 = np.sum(np.sqrt(moments[1] ** 2 + moments[2] ** 2))
+        err = e1 / e2
+        print('error is = ', err)
+    u0 = moments[1]
+    v0 = moments[2]
+
+print("total time taken: ", time.time()-start)
+
+# to plot contours of velocity
+x = np.linspace(0, 1, nx)
+y = np.linspace(0, 1, ny)
+X, Y = np.meshgrid(x, y)
+VelNorm = np.sqrt(moments[1, 1:-1, 1:-1] * moments[1, 1:-1, 1:-1] 
+                  + moments[2, 1:-1, 1:-1] * moments[2, 1:-1, 1:-1])
+
+plt.figure()
+# cp = plt.contourf(X, Y, np.transpose(VelNorm), 25, cmap=plt.cm.viridis)
+cp = plt.contour(X, Y, np.transpose(VelNorm), 15)
+cbar=plt.colorbar(cp)
+cbar.set_label('Velocity', fontsize=12)
+plt.xlabel('x', fontsize=12)
+plt.ylabel('y', fontsize=12)
+plt.xticks(fontsize=12)
+plt.yticks(fontsize=12)
+plt.savefig('velocity_cont1.png', bbox_inches='tight')
+

sim_params.py

@@ -0,0 +1,47 @@
+from d2q9 import D2Q9
+
+
+class SimParams(D2Q9):
+    def __init__(self, mach, reynolds_num, ref_len, res):
+        self.mach_num = mach
+        self.reynolds = reynolds_num
+        self.ref_len = ref_len
+        self.res = res
+        self.dt = self.res
+        D2Q9.__init__(self)
+
+    @property
+    def ref_vel(self):
+        """
+        Compute reference velocity
+        """
+        return self.mach_num * self.a
+
+    @property
+    def knudsen(self):
+        """
+        Compute Knudsen number
+        """
+        return self.mach_num / self.reynolds
+
+    @property
+    def nu(self):
+        """
+        Compute viscosity from reference length, velocity and Reynolds number
+        """
+        return self.ref_len * self.ref_vel / self.reynolds
+
+    @property
+    def tau0(self):
+        return self.nu / self.T0
+
+    @property
+    def tau_ndim(self):
+        return self.tau0 / self.dt
+
+    @property
+    def beta(self):
+        """
+        Relaxation 
+        """
+        return 1.0 / (2.0 * self.tau_ndim + 1.0)