Shakers

libd3q15/.gitignore

@@ -1 +1,2 @@
 eigenvectors.h
+config.mk

libd3q15/Makefile

@@ -15,6 +15,9 @@ libobjs = d3q15.o checks.o noise.o \
 
 headers = d3q15.h Lattice.h noise.h eigenvectors.h
 
+localconfig = $(wildcard config.mk)
+include $(localconfig)
+
 $(LIB) : $(libobjs)
 	$(AR) $@ $^
 	ranlib $@

libd3q15/d3q15.c

@@ -86,11 +86,11 @@ void d3q15_iterate(Lattice *lat, int n_steps) {
     collide(lat);
 
     /* Update boundary */
-    (*lat->bc_func)(lat);
+    (*lat->bc_func)(lat);     
 
     /* Propagate */
     propagate(lat);
-    
+
     lat->time_step++;
   }
 }
@@ -268,7 +268,7 @@ void collide (Lattice *lat) {
 	/* relax the trace */
 	TrS -= omega_b*(TrS - site.rho[0]*usq);
 	/* Add forcing part to trace */
-	TrS += 2.*omega_b*tau_b * uDOTf;
+	TrS += 2.*omega_b*tau_b * uDOTf / site.rho[0];
 
 	/* and the traceless part */
 	for (a=0; a<DQ_d; a++) {
@@ -277,7 +277,7 @@ void collide (Lattice *lat) {
 				site.rho[0]*(site.u[a]*site.u[b] -usq*delta[a][b]));
 
 	    /* including traceless force */
-	    S[a][b] += 2.*omega_s*tau_s * (site.u[a]*site.force[b] + site.force[a]*site.u[b] - 2. * uDOTf * delta[a][b]);
+  	    S[a][b] += omega_s*tau_s * (site.u[a]*site.force[b] + site.force[a]*site.u[b] - 2. * uDOTf * delta[a][b]) / site.rho[0];
 	  }
 	  /* add the trace back on */
 	  S[a][a] += TrS / DQ_d;

python/.gitignore

@@ -1,3 +1,3 @@
 *.pyc
 build
-
+setup.cfg

python/d3q15/Lattice.i

@@ -315,10 +315,21 @@ EXC_CHECK(force_set)
    */
   void initFromHydroVars() {
     int i, j, k;
+    int a;
+    double u[DQ_d], force[DQ_d];
+
     for (i=1; i<=$self->nx; i++)
       for (j=1; j<=$self->ny; j++)
-	for (k=1; k<=$self->nz; k++)
-	  calc_equil(DQ_rho_get($self,i,j,k), &DQ_u_get($self, i,j,k,0), &DQ_f_get($self, i,j,k,0));
+	for (k=1; k<=$self->nz; k++) {
+	  for (a=0; a<DQ_d; a++) {
+            u[a] = DQ_u_get($self, i,j,k, a);
+            force[a] = DQ_force_get($self, i,j,k, a);
+            u[a] += 0.5*force[a];
+/*          if(u[a] != 0.0) printf("%5i %5i %5i %5i %15.10e \n", i,j,k, a, u[a]); */
+          }  
+/*	  calc_equil(DQ_rho_get($self,i,j,k), &DQ_u_get($self, i,j,k,0), &DQ_f_get($self, i,j,k,0)); */
+	  calc_equil(DQ_rho_get($self,i,j,k), u, &DQ_f_get($self, i,j,k,0));
+        }
   }
 
   /* pseduo method wrapper */

python/d3q15/swimmers.py

@@ -65,15 +65,24 @@ def __init__(self, L, **kwargs):
         except KeyError:
             randSeed = 0
 
-        for flag in ['walls', 'tumbles', 'sinks', 'tracks', 'forced']:
+        # Set default options
+        self.__setattr__('swims', True)
+        self.__setattr__('advects', True)
+        self.__setattr__('rotates', True)
+        self.__setattr__('tumbles', True)
+        self.__setattr__('tracks', True)
+        self.__setattr__('walls', False)
+        self.__setattr__('sinks', False)
+        self.__setattr__('forced', False)
+
+        for flag in ['swims', 'advects', 'rotates', 'tumbles', 'tracks', 'walls', 'sinks', 'forced']:
             try:
                 if kwargs[flag]:
                     self.__setattr__(flag, True)
                 else:
                     self.__setattr__(flag, False)
             except KeyError:
-                self.__setattr__(flag, False)
-            continue
+                continue
 
         if self.tumbles:
             try: 
@@ -89,7 +98,7 @@ def __init__(self, L, **kwargs):
 
         self.eta = L.tau_s / 3.
         self.num = len(r_list)
-        
+
         # set up fields in data array
         # we only create a table if one doesn't exist already
         # (we might be called by a subclass that created some already)
@@ -167,9 +176,10 @@ def move(self, lattice):
         where v(R) is the interpolated velocity at the position of the
         swimmer, a is the radius and n_ is the orientation.
         """
-        v = self._interp(lattice, self.r)
-
-        rDot = self.P * self.n 
+        rDot = N.zeros([self.num,3])   
+        if self.swims:
+            rDot += self.P * self.n 
+            pass
         if self.sinks:
             rDot += self.F 
             pass
@@ -178,18 +188,24 @@ def move(self, lattice):
             pass
 
         h = self.hydroRadius
-
         rDot *= (1./self.a - 1./h) / (6. * N.pi * self.eta)
-        rDot += v
+
+        v = self._interp(lattice, self.r)
+        if self.advects:
+            rDot += v
+            pass
 
         # rPlus = self.r
-        rMinus= self.r - self.n * self.l
-        # from above, v = v(rPlus)
-        nDot = v - self._interp(lattice, rMinus)
-        # now nDot = v(rPlus) - v(rMinus)
-        # so divide by l to get the rate
-        nDot /= self.l
-
+        nDot = N.zeros([self.num,3])   
+        if self.rotates: 
+            rMinus = self.r - self.n * self.l
+            # from above, v = v(rPlus)
+            nDot += v - self._interp(lattice, rMinus)
+            # now nDot = v(rPlus) - v(rMinus)
+            # so divide by l to get the rate
+            nDot /= self.l
+            pass
+
         self.applyMove(lattice, rDot)
         self.applyTurn(lattice, nDot)
 
@@ -269,7 +285,7 @@ def applyMove(self, lattice, rdots):
             if self.walls:
                 # make sure the unwrapped coords do respect the walls
                 # by copying in the actual z
-                self.s[2] = self.r[2]
+                self.s[:,2] = self.r[:,2]
         return
 
     def applyTurn(self, lattice, nDots):

python/setup.py

@@ -5,6 +5,7 @@
 import distutils.sysconfig
 from Cython.Build import cythonize
 import numpy
+import numpy.distutils.intelccompiler
 
 numpyIncludeDir = numpy.get_include()