Lots of work on website and examples

examples/MD_CFD/cfd_cpl.py

@@ -15,7 +15,7 @@
 # Parameters of the cpu topology (cartesian grid)
 npxyz = np.array([1, 1, 1], order='F', dtype=np.int32)
 NProcs = np.product(npxyz)
-xyzL = np.array([10.0, 16.0, 1.0], order='F', dtype=np.float64)
+xyzL = np.array([6.70820393, 17.88854382, 1.0], order='F', dtype=np.float64)
 xyz_orig = np.array([0.0, 0.0, 0.0], order='F', dtype=np.float64)
 ncxyz = np.array([8, 8, 1], order='F', dtype=np.int32)
 
@@ -42,15 +42,17 @@
 A_send = np.zeros((2, ncxl, ncyl, nczl), order='F', dtype=np.float64)
 
 #Set CFD simulation object
-t0 = 0; tf = 30.; Nsteps = 2500
-time = np.linspace(t0, tf, Nsteps)
-dt = 0.2
+md_cfd_dt_ratio = 25
+dt = 0.125; Nsteps = 100000/md_cfd_dt_ratio; tf = Nsteps*dt
+time = np.arange(-dt,tf,dt)
 uwall = 1.
 cfd = CFD(nu=0.575, dt=dt, 
           xsize = ncxyz[0], ysize = ncxyz[1]+2,
           xmin = xyz_orig[0], xmax = xyzL[0],
           ymin = xyz_orig[1], ymax = xyzL[1])
 
+
+
 #Main Run
 for n,t in enumerate(time):
 
@@ -60,23 +62,21 @@
     # Call to CPL-LIBRARY goes here to
     # send u_CFD in constraint region
     #===============================================
-    A_send[0,:,0,0] = cfd.u[:,-3]
+    A_send[0,:,0,0] = cfd.u[:,2]
     CPL.send(A_send, limits_MD_BC)
 
     #===============================================
     # Call to CPL-LIBRARY goes here to
     # recieve u_MD to set bottom boundary
     #===============================================
     umd, ierr = CPL.recv(A_recv, limits_CFD_BC)
-    bottomwall = np.mean(umd)
-
-    #Update CFD
+    bottomwall = umd[0,:,0,0]
     cfd.set_bc(topwall=uwall, bottomwall=bottomwall)
-    cfd.update_time()
 
     #plot
     cfd.plot()
 
-
+    #Update CFD timestep
+    cfd.update_time()
 
 

examples/MD_CFD/cfd_oo.py

@@ -61,10 +61,11 @@ def __init__(self, dt, nu = 1.,
             self.ax = fig.axes
 
     def set_bc(self, topwall=1., bottomwall=0.):
+        #Periodic boundaries
+        self.u[-1,1:-1] = self.u[1,1:-1]; self.u[0,1:-1] = self.u[-2,1:-1]
         #Enforce boundary conditions
         self.u[:,0] = bottomwall; self.u[:,-1] = topwall
-        #Periodic boundaries
-        self.u[-1,:] = self.u[1,:]; self.u[0,:] = self.u[-2,:]
+
 
     def update_time(self):
         # Save previous value
@@ -84,25 +85,35 @@ def plot(self, ax=None):
         if ax == None:
             ax=self.ax
 
-        sm = ax.pcolormesh(self.X,self.Y,self.u.T,vmin=0.,vmax=1.,alpha=0.5)
-        draw_grid(ax, nx=self.x.size-1,ny=self.y.size-1, nz=1,
+        sm = ax.imshow(self.u.T,aspect='auto',origin='lower',
+                       extent=[self.xmin,self.xmax,
+                               self.ymin,self.ymax],
+                       interpolation="none",vmin=-1.,vmax=1.,
+                       alpha=0.5, cmap=plt.cm.RdYlBu_r)
+
+#        sm = ax.pcolormesh(self.X,self.Y,self.u.T,vmin=-1.,vmax=1.,alpha=0.5,
+#                          cmap=plt.cm.RdYlBu_r)
+        draw_grid(ax, nx=self.x.size,ny=self.y.size, nz=1,
                       xmin=self.xmin,xmax=self.xmax,
                       ymin=self.ymin,ymax=self.ymax)
 
         #Plot velocity profile offset to the left
-        axisloc = self.xmax
-        ax.arrow(axisloc,0., 0., 1.,  width=0.0015, color="k", 
-                 clip_on=False, head_width=0.012, head_length=0.012)
-        ax.arrow(axisloc,0., 1., 0.,  width=0.0015, color="k", 
-                 clip_on=False, head_width=0.012, head_length=0.012)
-        ax.plot(np.mean(self.u,0)+axisloc,self.y,'g-')
+        axisloc = self.xmax+1.
+
+        ax.arrow(axisloc, 0., self.ymin, self.ymax,  width=0.0015, color="k", 
+                 clip_on=False, head_width=0.12, head_length=0.12)
+        mid = .5*(self.ymin+self.ymax)
+        ax.arrow(axisloc-1., mid, 2.0, 0., width=0.0015, color="k", 
+                 clip_on=False, head_width=0.12, head_length=0.12)
+        yp = np.linspace(self.ymin+.5*self.dy, self.ymax - 0.5*self.dy, self.y.size)
+        ax.plot(np.mean(self.u,0)+axisloc,yp,'g-x')
         #ax.set_xlim((0.,2.))
 
         if self.first_time:
             plt.colorbar(sm)
             self.first_time=False
 
-        plt.pause(0.01)
+        plt.pause(0.001)
         ax.cla()
 
 if __name__ == "__main__":

examples/MD_CFD/md_cpl.py

@@ -15,7 +15,7 @@
 # Parameters of the cpu topology (cartesian grid)
 npxyz = np.array([1, 1, 1], order='F', dtype=np.int32)
 NProcs = np.product(npxyz)
-xyzL = np.array([10.0, 16.0, 1.0], order='F', dtype=np.float64)
+xyzL = np.array([6.70820393, 17.88854382, 1.0], order='F', dtype=np.float64)
 xyz_orig = np.array([0.0, 0.0, 0.0], order='F', dtype=np.float64)
 
 if (nprocs_realm != NProcs):
@@ -42,15 +42,15 @@
 A_recv = np.zeros((2, ncxl, ncyl, nczl), order='F', dtype=np.float64)
 
 # Setup MD simulation object
-t0 = 0; tf = 30.; Nsteps = 10000
-time = np.linspace(t0, tf, Nsteps)
-dt = 0.005
+md_cfd_dt_ratio = 25
+dt = 0.005; Nsteps = 100000; tf = Nsteps*dt
+time = np.arange(0.,tf,dt)
 md = MD(dt=dt, wallwidth=[2.,0.], wallslide=[-1.,0.])
 
 #Main run
 for n,t in enumerate(time):
 
-    print("MD time = ", n,t)
+    print("MD time = ", md.tstep, md.time)
 
     # Calculate force
     md.force()
@@ -61,24 +61,30 @@
     # and force is applied
     # F = (1/tau)*(u_CFD - u_MD)
     #=======================================================
-    A_recv, ierr = CPL.recv(A_recv, limits_MD_BC)
-    u_CFD = A_recv[0,:,:,:]
-    print(A_recv.shape, u_CFD.shape)
-    md.constraint_force(u_CFD, 8)
+    if n%md_cfd_dt_ratio == 0:        
+        A_recv, ierr = CPL.recv(A_recv, limits_MD_BC)
+        u_CFD = A_recv[:,:,:,0]
+
+    #Cell 7 is constrained
+    md.constraint_force(u_CFD, 7)
 
     # Calculate velocity
     md.verlet()
 
+    #Plot
+    if n%md_cfd_dt_ratio == 0:        
+        md.plot()
+
     #=======================================================
     #Call to CPL-LIBRARY to send u_MD at boundary
     #=======================================================
-    u = md.get_velfield([ncx,ncy])
-    A_send[0,:,0,0] = u[:,6,0] #Cell 6 is 
-    A_send[1,:,0,0] = u[:,6,1]
-    CPL.send(A_send, limits_CFD_BC)
+    if n%md_cfd_dt_ratio == 0:        
+        u = md.get_velfield([ncx,ncy])
+        #Cell 5 is sent
+        A_send[0,:,0,0] = u[0,:,5]
+        A_send[1,:,0,0] = u[1,:,5]
+        CPL.send(A_send, limits_CFD_BC)
+
 
-    #Plot
-    md.plot()
-