diff --git a/pystokes/__init__.py b/pystokes/__init__.py index 263c31d..eb60d0e 100644 --- a/pystokes/__init__.py +++ b/pystokes/__init__.py @@ -1,5 +1,6 @@ import pystokes.interface import pystokes.periodic +import pystokes.nearest_neighbors import pystokes.twoWalls import pystokes.unbounded import pystokes.utils diff --git a/pystokes/forceFields.pxd b/pystokes/forceFields.pxd index 5f82863..38cb4cd 100644 --- a/pystokes/forceFields.pxd +++ b/pystokes/forceFields.pxd @@ -60,6 +60,8 @@ cdef class Forces: cpdef membraneSurface(self, int Nmx, int Nmy, double [:] F, double [:] r, double bondLength, double springModulus, double bendModulus ) + cpdef Cosserat(self, double [:] F, double [:] r, double [:] e1, double [:] e2, double [:] e3, double Lambda, double d) + @@ -73,3 +75,5 @@ cdef class Torques: cpdef bottomHeaviness(self, double [:] T, double [:] p, double bh=?) cpdef magnetic(self, double[:] T, double [:] p, double m0, double Bx, double By, double Bz) + + cpdef Cosserat(self, double [:] T, double [:] r, double [:] e1, double [:] e2, double [:] e3, double Lambda, double mu, double d) diff --git a/pystokes/forceFields.pyx b/pystokes/forceFields.pyx index 9357708..df3908f 100644 --- a/pystokes/forceFields.pyx +++ b/pystokes/forceFields.pyx @@ -984,7 +984,62 @@ cdef class Forces: return + cpdef Cosserat(self, double [:] F, double [:] r, double [:] e1, double [:] e2, double [:] e3, double Lambda, double d): + """ + The force for the Cosserat solids + ... + + Parameter + ---------- + r: np.array + An array of positions + An array of size 3*N, + F: np.array + An array of forces + An array of size 3*N, + e1, e2, e3: np.array + Arrays of Orientations + Arrays of size 3*N + Lambda: float + Strength of the force + d: float + particle distance + """ + cdef int N = self.N, i, k, xx = 2*N, ip, im + cdef double L = N*d + cdef double dxp, dxm, dyp, dym, dzp, dzm + + for i in prange(N,nogil=True): + ip = (i + 1) % N + im = (i - 1 + N) % N + + dxp = r[ip] - r[i] + dxm = r[i] - r[im] + dyp = r[N + ip] - r[N + i] + dym = r[N + i] - r[N + im] + dzp = r[xx + ip] - r[xx + i] + dzm = r[xx + i] - r[xx + im] + + if i == 0: + dxm = r[i] - r[im] + L + elif i == N - 1: + dxp = r[ip] - r[i] + L + + for k in range(3): + F[i+k*N] += 0.5 * Lambda * ((e1[i] * dxp + e1[N+i] * dyp + e1[xx+i] * dzp - d) * e1[i+k*N] \ + + (e2[i] * dxp + e2[N+i] * dyp + e2[xx+i] * dzp) * e2[i+k*N] \ + + (e3[i] * dxp + e3[N+i] * dyp + e3[xx+i] * dzp) * e3[i+k*N]) + F[i+k*N] -= 0.5 * Lambda * ((e1[i] * dxm + e1[N+i] * dym + e1[xx+i] * dzm - d) * e1[i+k*N] \ + + (e2[i] * dxm + e2[N+i] * dym + e2[xx+i] * dzm) * e2[i+k*N] \ + + (e3[i] * dxm + e3[N+i] * dym + e3[xx+i] * dzm) * e3[i+k*N]) + F[i+k*N] += 0.5 * Lambda * ((e1[ip] * dxp + e1[N+ip] * dyp + e1[xx+ip] * dzp - d) * e1[ip+k*N] \ + + (e2[ip] * dxp + e2[N+ip] * dyp + e2[xx+ip] * dzp) * e2[ip+k*N] \ + + (e3[ip] * dxp + e3[N+ip] * dyp + e3[xx+ip] * dzp) * e3[ip+k*N]) + F[i+k*N] -= 0.5 * Lambda * ((e1[im] * dxm + e1[N+im] * dym + e1[xx+im] * dzm - d) * e1[im+k*N] \ + + (e2[im] * dxm + e2[N+im] * dym + e2[xx+im] * dzm) * e2[im+k*N] \ + + (e3[im] * dxm + e3[N+im] * dym + e3[xx+im] * dzm) * e3[im+k*N]) + return @cython.boundscheck(False) @cython.cdivision(True) @@ -1067,3 +1122,94 @@ cdef class Torques: T[i+Z] += m0*(p[i]*By -p[i+N]*Bx) return + + cpdef Cosserat(self, double [:] T, double [:] r, double [:] e1, double [:] e2, double [:] e3, double Lambda, double mu, double d): + """ + The torque for the Cosserat solids + + ... + + Parameter + ---------- + r: np.array + An array of positions + An array of size 3*N, + T: np.array + An array of torques + An array of size 3*N, + e1, e2, e3: np.array + Arrays of Orientations + Arrays of size 3*N + Lambda: float + Strength of the force + mu: float + Strength of the torque + d: float + particle distance + """ + cdef int N = self.N, i, j, k, xx = 2*N, ip, im + cdef double L = N*d + cdef double drp[3], drm[3], e[3][3], ep[3][3], em[3][3] + cdef double e_dot_r, e_dot_e + cdef double cross[3] + + for i in prange(N,nogil=True): + ip = (i + 1) % N + im = (i - 1 + N) % N + + drp[0] = r[ip] - r[i] + drp[1] = r[N+ip] - r[N+i] + drp[2] = r[xx + ip] - r[xx + i] + + drm[0] = r[i] - r[im] + drm[1] = r[N + i] - r[N + im] + drm[2] = r[xx + i] - r[xx + im] + + if i == 0: + drm[0] += L + elif i == N - 1: + drp[0] += L + + for k in range(3): + e[0][k] = e1[k*N+i] + e[1][k] = e2[k*N+i] + e[2][k] = e3[k*N+i] + + ep[0][k] = e1[k*N+ip] + ep[1][k] = e2[k*N+ip] + ep[2][k] = e3[k*N+ip] + + em[0][k] = e1[k*N+im] + em[1][k] = e2[k*N+im] + em[2][k] = e3[k*N+im] + + for j in range(3): + e_dot_r = e[j][0] * drp[0] + e[j][1] * drp[1] + e[j][2] * drp[2] + cross[0] = e[j][1] * drp[2] - e[j][2] * drp[1] + cross[1] = e[j][2] * drp[0] - e[j][0] * drp[2] + cross[2] = e[j][0] * drp[1] - e[j][1] * drp[0] + for k in range(3): + T[i+k*N] -= 0.5 * Lambda * (e_dot_r - (d if j == 0 else 0)) * cross[k] + + e_dot_r = e[j][0] * drm[0] + e[j][1] * drm[1] + e[j][2] * drm[2] + cross[0] = e[j][1] * drm[2] - e[j][2] * drm[1] + cross[1] = e[j][2] * drm[0] - e[j][0] * drm[2] + cross[2] = e[j][0] * drm[1] - e[j][1] * drm[0] + for k in range(3): + T[i+k*N] -= 0.5 * Lambda * (e_dot_r - (d if j == 0 else 0)) * cross[k] + + for j in range(3): + e_dot_e = e[j][0] * ep[(j+1)%3][0] + e[j][1] * ep[(j+1)%3][1] + e[j][2] * ep[(j+1)%3][2] + cross[0] = e[j][1] * ep[(j+1)%3][2] - e[j][2] * ep[(j+1)%3][1] + cross[1] = e[j][2] * ep[(j+1)%3][0] - e[j][0] * ep[(j+1)%3][2] + cross[2] = e[j][0] * ep[(j+1)%3][1] - e[j][1] * ep[(j+1)%3][0] + for k in range(3): + T[i+k*N] -= 0.5 * mu * (e_dot_e) * cross[k] + + e_dot_e = e[j][0] * em[(j+2)%3][0] + e[j][1] * em[(j+2)%3][1] + e[j][2] * em[(j+2)%3][2] + cross[0] = e[j][1] * em[(j+2)%3][2] - e[j][2] * em[(j+2)%3][1] + cross[1] = e[j][2] * em[(j+2)%3][0] - e[j][0] * em[(j+2)%3][2] + cross[2] = e[j][0] * em[(j+2)%3][1] - e[j][1] * em[(j+2)%3][0] + for k in range(3): + T[i+k*N] -= 0.5 * mu * (e_dot_e) * cross[k] + return \ No newline at end of file diff --git a/pystokes/nearest_neighbors.pxd b/pystokes/nearest_neighbors.pxd new file mode 100644 index 0000000..8dd2a0b --- /dev/null +++ b/pystokes/nearest_neighbors.pxd @@ -0,0 +1,41 @@ +cimport cython +from libc.math cimport sqrt +from cython.parallel import prange +import numpy as np +cimport numpy as np +cdef double PI = 3.14159265359 + +@cython.wraparound(False) +@cython.boundscheck(False) +@cython.cdivision(True) +@cython.nonecheck(False) +cdef class Rbm: + cdef double a, eta, L, mu, muv, mur + cdef int N + cdef readonly np.ndarray Mobility + + cpdef mobilityTT(self, double [:] v, double [:] r, double [:] F) + + cpdef mobilityTR(self, double [:] v, double [:] r, double [:] T) + + cpdef propulsionT2s(self, double [:] v, double [:] r, double [:] S) + + cpdef propulsionT3t(self, double [:] v, double [:] r, double [:] D) + + cpdef propulsionT3a(self, double [:] v, double [:] r, double [:] V) + + cpdef propulsionT4a(self, double [:] v, double [:] r, double [:] M) + + + ## Angular velocities + + + cpdef mobilityRT(self, double [:] o, double [:] r, double [:] F) + + cpdef mobilityRR( self, double [:] o, double [:] r, double [:] T) + + cpdef propulsionR2s(self, double [:] o, double [:] r, double [:] S) + + cpdef propulsionR3a(self, double [:] v, double [:] r, double [:] V) + + cpdef propulsionR4a( self, double [:] o, double [:] r, double [:] M) \ No newline at end of file diff --git a/pystokes/nearest_neighbors.pyx b/pystokes/nearest_neighbors.pyx new file mode 100644 index 0000000..598e44a --- /dev/null +++ b/pystokes/nearest_neighbors.pyx @@ -0,0 +1,659 @@ +cimport cython +from libc.math cimport sqrt +from cython.parallel import prange +cdef double PI = 3.14159265359 +import numpy as np +cimport numpy as np + + +@cython.wraparound(False) +@cython.boundscheck(False) +@cython.cdivision(True) +@cython.nonecheck(False) +cdef class Rbm: + """ + Rigid body motion (RBM) - velocity and angular velocity + + Methods in this class update velocities or angular velocities + using the inputs of - arrays of positions, velocity or angular velocity, + along with an array of forces or torques or a slip mode + + The array of velocity or angular velocities is then update by each method. + + ... + + ---------- + radius: float + Radius of the particles (a). + particles: int + Number of particles (N) + viscosity: float + Viscosity of the fluid (eta) + boxSize: float + Length of the box which is reperated periodicly in 3D + """ + + def __init__(self, radius=1, particles=1, viscosity=1.0, boxSize=10.0): + self.a = radius + self.N = particles + self.eta = viscosity + self.L = boxSize + self.mu = 1.0/(6*PI*self.eta*self.a) + self.muv = 1.0/(8*PI*self.eta) + self.mur = 1.0/(8*PI*self.eta*self.a**3) + + self.Mobility = np.zeros( (3*self.N, 3*self.N), dtype=np.float64) + + + cpdef mobilityTT(self, double [:] v, double [:] r, double [:] F): + """ + Compute velocity due to body forces using :math:`v=\mu^{TT}\cdot F` + ... + + Parameters + ---------- + v: np.array + An array of velocities + An array of size 3*N, + r: np.array + An array of positions + An array of size 3*N, + F: np.array + An array of forces + An array of size 3*N, + """ + + + cdef int N = self.N, i, j, xx=2*N + cdef double dx, dy, dz, idr, idr2, vx, vy, vz, vv1, vv2, aa = (2.0*self.a*self.a)/3.0, L = self.L + cdef double mu=self.mu, muv=self.muv + cdef int neighbors[2] + + for i in prange(N, nogil=True): + vx=0; vy=0; vz=0; + neighbors[0] = i - 1 + neighbors[1] = i + 1 + for j in neighbors: + if i == 0 and j == i - 1: + j = N - 1 + dx = r[i] - (r[j] - L) + elif i == N - 1 and j == i + 1: + j = 0 + dx = r[i] - (r[j] + L) + else: + dx = r[i] - r[j] + dy = r[i+N] - r[j+N] + dz = r[i+xx] - r[j+xx] + idr = 1.0/sqrt( dx*dx + dy*dy + dz*dz ) + idr2 = idr*idr + + vv1 = (1+aa*idr2)*idr + vv2 = (1-3*aa*idr2)*( F[j]*dx + F[j+N]*dy + F[j+xx]*dz )*idr2*idr + vx += vv1*F[j] + vv2*dx + vy += vv1*F[j+N] + vv2*dy + vz += vv1*F[j+xx] + vv2*dz + + v[i] += mu*F[i] + muv*vx + v[i+N] += mu*F[i+N] + muv*vy + v[i+xx] += mu*F[i+xx] + muv*vz + return + + + cpdef mobilityTR(self, double [:] v, double [:] r, double [:] T): + """ + Compute velocity due to body torque using :math:`v=\mu^{TR}\cdot T` + ... + + Parameters + ---------- + v: np.array + An array of velocities + An array of size 3*N, + r: np.array + An array of positions + An array of size 3*N, + T: np.array + An array of torques + An array of size 3*N, + """ + + + cdef int N = self.N, i, j, xx=2*N + cdef double dx, dy, dz, idr, idr3, vx, vy, vz, L = self.L + cdef double muv=self.muv + cdef int neighbors[2] + + for i in prange(N, nogil=True): + vx=0; vy=0; vz=0; + neighbors[0] = i - 1 + neighbors[1] = i + 1 + for j in neighbors: + if i == 0 and j == i - 1: + j = N - 1 + dx = r[i] - (r[j] - L) + elif i == N - 1 and j == i + 1: + j = 0 + dx = r[i] - (r[j] + L) + else: + dx = r[i] - r[j] + dy = r[i+N] - r[j+N] + dz = r[i+xx] - r[j+xx] + idr = 1.0/sqrt( dx*dx + dy*dy + dz*dz ) + idr3 = idr*idr*idr + vx += -(dy*T[j+xx] -T[j+N]*dz )*idr3 + vy += -(dz*T[j] -T[j+xx]*dx )*idr3 + vz += -(dx*T[j+N] -T[j] *dy )*idr3 + + v[i] += muv*vx + v[i+N] += muv*vy + v[i+xx] += muv*vz + return + + + cpdef propulsionT2s(self, double [:] v, double [:] r, double [:] S): + """ + Compute velocity due to 2s mode of the slip :math:`v=\pi^{T,2s}\cdot S` + ... + + Parameters + ---------- + v: np.array + An array of velocities + An array of size 3*N, + r: np.array + An array of positions + An array of size 3*N, + S: np.array + An array of 2s mode of the slip + An array of size 5*N, + """ + + cdef int N = self.N, i, j, xx=2*N, xx1=3*N, xx2=4*N + cdef double dx, dy, dz, dr, idr, idr3, L = self.L + cdef double aa=(self.a*self.a*8.0)/3.0, vv1, vv2, aidr2 + cdef double vx, vy, vz, + cdef double sxx, sxy, sxz, syz, syy, srr, srx, sry, srz, mus = (28.0*self.a*self.a)/24 + cdef int neighbors[2] + + for i in prange(N, nogil=True): + vx=0; vy=0; vz=0; + neighbors[0] = i - 1 + neighbors[1] = i + 1 + for j in neighbors: + if i == 0 and j == i - 1: + j = N - 1 + dx = r[i] - (r[j] - L) + elif i == N - 1 and j == i + 1: + j = 0 + dx = r[i] - (r[j] + L) + else: + dx = r[i] - r[j] + sxx = S[j] + syy = S[j+N] + sxy = S[j+xx] + sxz = S[j+xx1] + syz = S[j+xx2] + dy = r[i+N] - r[j+N] + dz = r[i+xx] - r[j+xx] + idr = 1.0/sqrt( dx*dx + dy*dy + dz*dz ) + idr3 = idr*idr*idr + aidr2 = aa*idr*idr + + srr = (sxx*(dx*dx-dz*dz) + syy*(dy*dy-dz*dz) + 2*sxy*dx*dy + 2*sxz*dx*dz + 2*syz*dy*dz)*idr*idr + srx = sxx*dx + sxy*dy + sxz*dz + sry = sxy*dx + syy*dy + syz*dz + srz = sxz*dx + syz*dy - (sxx+syy)*dz + + vv1 = 3*(1-aidr2)*srr*idr3 + vv2 = 1.2*aidr2*idr3 + vx += vv1*dx + vv2*srx + vy += vv1*dy + vv2*sry + vz += vv1*dz + vv2*srz + + v[i] += vx*mus + v[i+N]+= vy*mus + v[i+xx]+= vz*mus + + return + + cpdef propulsionT3t(self, double [:] v, double [:] r, double [:] D): + """ + Compute velocity due to 3t mode of the slip :math:`v=\pi^{T,3t}\cdot D` + ... + + Parameters + ---------- + v: np.array + An array of velocities + An array of size 3*N, + r: np.array + An array of positions + An array of size 3*N, + D: np.array + An array of 3t mode of the slip + An array of size 3*N, + """ + + cdef int N = self.N, i, j, xx=2*N + cdef double dx, dy, dz, idr, idr3, Ddotidr, vx, vy, vz, mud = 3.0*self.a*self.a*self.a/5, mud1 = -1.0*(self.a**3)/5 + cdef double L = self.L + cdef int neighbors[2] + + for i in prange(N, nogil=True): + vx=0; vy=0; vz=0; + neighbors[0] = i - 1 + neighbors[1] = i + 1 + for j in neighbors: + if i == 0 and j == i - 1: + j = N - 1 + dx = r[i] - (r[j] - L) + elif i == N - 1 and j == i + 1: + j = 0 + dx = r[i] - (r[j] + L) + else: + dx = r[i] - r[j] + + dy = r[i+N] - r[j+N] + dz = r[i+xx] - r[j+xx] + idr = 1.0/sqrt( dx*dx + dy*dy + dz*dz ) + idr3 = idr*idr*idr + Ddotidr = (D[j]*dx + D[j+N]*dy + D[j+xx]*dz)*idr*idr + + vx += (D[j] - 3.0*Ddotidr*dx )*idr3 + vy += (D[j+N] - 3.0*Ddotidr*dy )*idr3 + vz += (D[j+xx] - 3.0*Ddotidr*dz )*idr3 + + v[i] += mud1*vx + v[i+N] += mud1*vy + v[i+xx]+= mud1*vz + return + + cpdef propulsionT3a(self, double [:] v, double [:] r, double [:] V): + """ + Compute velocity due to 3a mode of the slip :math:`v=\pi^{T,3a}\cdot V` + ... + + Parameters + ---------- + v: np.array + An array of velocities + An array of size 3*N, + r: np.array + An array of positions + An array of size 3*N, + V: np.array + An array of 3a mode of the slip + An array of size 5*N, + """ + + cdef int N = self.N, i, j + cdef double dx, dy, dz, idr, idr5, vxx, vyy, vxy, vxz, vyz, vrx, vry, vrz + cdef double mud = 13.0*self.a*self.a*self.a/12, L = self.L + cdef int neighbors[2] + + for i in prange(N, nogil=True): + neighbors[0] = i - 1 + neighbors[1] = i + 1 + for j in neighbors: + if i == 0 and j == i - 1: + j = N - 1 + dx = r[i] - (r[j] - L) + elif i == N - 1 and j == i + 1: + j = 0 + dx = r[i] - (r[j] + L) + else: + dx = r[i] - r[j] + vxx = V[j] + vyy = V[j+N] + vxy = V[j+2*N] + vxz = V[j+3*N] + vyz = V[j+4*N] + + dy = r[i+N] - r[j+N] + dz = r[i+2*N] - r[j+2*N] + idr = 1.0/sqrt( dx*dx + dy*dy + dz*dz) + idr5 = idr*idr*idr*idr*idr + vrx = vxx*dx + vxy*dy + vxz*dz + vry = vxy*dx + vyy*dy + vyz*dz + vrz = vxz*dx + vyz*dy - (vxx+vyy)*dz + + v[i] += mud * (vry * dz - vrz * dy) * idr5 + v[i+N] += mud * (vrz * dx - vrx * dz) * idr5 + v[i+2*N] += mud * (vrx * dy - vry * dx) * idr5 + return + + cpdef propulsionT4a(self, double [:] v, double [:] r, double [:] M): + """ + Compute velocity due to 4a mode of the slip :math:`v=\pi^{T,4a}\cdot M` + ... + + Parameters + ---------- + v: np.array + An array of velocities + An array of size 3*N, + r: np.array + An array of positions + An array of size 3*N, + M: np.array + An array of 4a mode of the slip + An array of size 7*N, + """ + + cdef int N = self.N, i, j + cdef double dx, dy, dz, idr, idr7 + cdef double mrrx, mrry, mrrz, mxxx, myyy, mxxy, mxxz, mxyy, mxyz, myyz + cdef double mud = -363/8 * self.a ** 4, L = self.L + cdef int neighbors[2] + + for i in prange(N, nogil=True): + neighbors[0] = i - 1 + neighbors[1] = i + 1 + for j in neighbors: + if i == 0 and j == i - 1: + j = N - 1 + dx = r[i] - (r[j] - L) + elif i == N - 1 and j == i + 1: + j = 0 + dx = r[i] - (r[j] + L) + else: + dx = r[i] - r[j] + + mxxx = M[j] + myyy = M[j+N] + mxxy = M[j+2*N] + mxxz = M[j+3*N] + mxyy = M[j+4*N] + mxyz = M[j+5*N] + myyz = M[j+6*N] + + dy = r[i+N] - r[j+N] + dz = r[i+2*N] - r[j+2*N] + idr = 1.0/sqrt( dx*dx + dy*dy + dz*dz ) + idr7 = idr*idr*idr*idr*idr*idr*idr + mrrx = mxxx*(dx*dx-dz*dz) + mxyy*(dy*dy-dz*dz) + 2*mxxy*dx*dy + 2*mxxz*dx*dz + 2*mxyz*dy*dz + mrry = mxxy*(dx*dx-dz*dz) + myyy*(dy*dy-dz*dz) + 2*mxyy*dx*dy + 2*mxyz*dx*dz + 2*myyz*dy*dz + mrrz = mxxz*(dx*dx-dz*dz) + myyz*(dy*dy-dz*dz) + 2*mxyz*dx*dy - 2*(mxxx+mxyy)*dx*dz - 2*(mxxy+myyy)*dy*dz + + v[i] -= mud*( dy*mrrz - dz*mrry )*idr7 + v[i+N] -= mud*( dz*mrrx - dx*mrrz )*idr7 + v[i+2*N] -= mud*( dx*mrry - dy*mrrx )*idr7 + return + + + ## Angular velocities + cpdef mobilityRT(self, double [:] o, double [:] r, double [:] F): + """ + Compute angular velocity due to body forces using :math:`o=\mu^{RT}\cdot F` + ... + + Parameters + ---------- + o: np.array + An array of angular velocities + An array of size 3*N, + r: np.array + An array of positions + An array of size 3*N, + F: np.array + An array of forces + An array of size 3*N, + """ + + cdef int N = self.N, i, j, xx=2*N + cdef double dx, dy, dz, idr, idr3, ox, oy, oz, muv=self.muv, L = self.L + cdef int neighbors[2] + + for i in prange(N, nogil=True): + ox=0; oy=0; oz=0; + neighbors[0] = i - 1 + neighbors[1] = i + 1 + for j in neighbors: + if i == 0 and j == i - 1: + j = N - 1 + dx = r[i] - (r[j] - L) + elif i == N - 1 and j == i + 1: + j = 0 + dx = r[i] - (r[j] + L) + else: + dx = r[i] - r[j] + dy = r[i+N] - r[j+N] + dz = r[i+xx] - r[j+xx] + idr = 1.0/sqrt( dx*dx + dy*dy + dz*dz ) + idr3 = idr*idr*idr + + ox += (F[j+N]*dz - F[j+xx]*dy )*idr3 + oy += (F[j+xx]*dx - F[j] *dz )*idr3 + oz += (F[j] *dy - F[j+N]*dx )*idr3 + o[i] += muv*ox + o[i+N] += muv*oy + o[i+xx] += muv*oz + return + + + cpdef mobilityRR( self, double [:] o, double [:] r, double [:] T): + """ + Compute angular velocity due to body torques using :math:`o=\mu^{RR}\cdot T` + ... + + Parameters + ---------- + o: np.array + An array of angular velocities + An array of size 3*N, + r: np.array + An array of positions + An array of size 3*N, + T: np.array + An array of forces + An array of size 3*N, + """ + + cdef int N = self.N, i, j, xx=2*N + cdef double dx, dy, dz, idr, idr3, Tdotidr, ox, oy, oz, mur=self.mur, muv=self.muv, L = self.L + cdef int neighbors[2] + + for i in prange(N, nogil=True): + ox=0; oy=0; oz=0; + neighbors[0] = i - 1 + neighbors[1] = i + 1 + for j in neighbors: + if i == 0 and j == i - 1: + j = N - 1 + dx = r[i] - (r[j] - L) + elif i == N - 1 and j == i + 1: + j = 0 + dx = r[i] - (r[j] + L) + else: + dx = r[i] - r[j] + dy = r[i+N] - r[j+N] + dz = r[i+xx] - r[j+xx] + idr = 1.0/sqrt( dx*dx + dy*dy + dz*dz ) + idr3 = idr*idr*idr + Tdotidr = ( T[j]*dx + T[j+N]*dy + T[j+xx]*dz )*idr*idr + + ox += ( T[j] - 3*Tdotidr*dx )*idr3 + oy += ( T[j+N] - 3*Tdotidr*dy )*idr3 + oz += ( T[j+xx] - 3*Tdotidr*dz )*idr3 + + o[i] += mur*T[i] - 0.5*muv*ox ##changed factor 0.5 here + o[i+N] += mur*T[i+N] - 0.5*muv*oy + o[i+xx] += mur*T[i+xx] - 0.5*muv*oz + return + + + cpdef propulsionR2s(self, double [:] o, double [:] r, double [:] S): + """ + Compute angular velocity due to 2s mode of the slip :math:`v=\pi^{R,2s}\cdot S` + ... + + Parameters + ---------- + o: np.array + An array of angular velocities + An array of size 3*N, + r: np.array + An array of positions + An array of size 3*N, + S: np.array + An array of 2s mode of the slip + An array of size 5*N, + """ + + cdef int N = self.N, i, j, xx=2*N + cdef double dx, dy, dz, idr, idr5, ox, oy, oz, L = self.L + cdef double sxx, sxy, sxz, syz, syy, srr, srx, sry, srz, mus = (28.0*self.a*self.a)/24 + cdef int neighbors[2] + + for i in prange(N, nogil=True): + ox=0; oy=0; oz=0; + neighbors[0] = i - 1 + neighbors[1] = i + 1 + for j in neighbors: + if i == 0 and j == i - 1: + j = N - 1 + dx = r[i] - (r[j] - L) + elif i == N - 1 and j == i + 1: + j = 0 + dx = r[i] - (r[j] + L) + else: + dx = r[i] - r[j] + sxx = S[j] + syy = S[j+N] + sxy = S[j+2*N] + sxz = S[j+3*N] + syz = S[j+4*N] + dy = r[i+N] - r[j+N] + dz = r[i+2*N] - r[j+2*N] + idr = 1.0/sqrt( dx*dx + dy*dy + dz*dz ) + idr5 = idr*idr*idr*idr*idr + srx = sxx*dx + sxy*dy + sxz*dz + sry = sxy*dx + syy*dy + syz*dz + srz = sxz*dx + syz*dy - (sxx+syy)*dz + + ox += 3*(sry*dz - srz*dy )*idr5 + oy += 3*(srz*dx - srx*dz )*idr5 + oz += 3*(srx*dy - sry*dx )*idr5 + + o[i] += ox*mus + o[i+N] += oy*mus + o[i+xx] += oz*mus + return + + cpdef propulsionR3a( self, double [:] o, double [:] r, double [:] V): + """ + Compute angular velocity due to 3a mode of the slip :math:`v=\pi^{R,3a}\cdot V` + ... + + Parameters + ---------- + o: np.array + An array of angular velocities + An array of size 3*N, + r: np.array + An array of positions + An array of size 3*N, + V: np.array + An array of 3a mode of the slip + An array of size 5*N, + """ + + cdef int N = self.N, i, j + cdef double dx, dy, dz, idr, idr2, idr5, vxx, vyy, vxy, vxz, vyz, vrr, vrx, vry, vrz + cdef double mud = 13.0 / 24.0 * self.a * self.a * self.a, L = self.L + cdef int neighbors[2] + + for i in prange(N, nogil=True): + neighbors[0] = i - 1 + neighbors[1] = i + 1 + for j in neighbors: + if i == 0 and j == i - 1: + j = N - 1 + dx = r[i] - (r[j] - L) + elif i == N - 1 and j == i + 1: + j = 0 + dx = r[i] - (r[j] + L) + else: + dx = r[i] - r[j] + + vxx = V[j] + vyy = V[j+N] + vxy = V[j+2*N] + vxz = V[j+3*N] + vyz = V[j+4*N] + + dy = r[i+N] - r[j+N] + dz = r[i+2*N] - r[j+2*N] + idr = 1.0/sqrt( dx*dx + dy*dy + dz*dz ) + idr5 = idr*idr*idr*idr*idr + vrr = (vxx*(dx*dx-dz*dz) + vyy*(dy*dy-dz*dz) + 2*vxy*dx*dy + 2*vxz*dx*dz + 2*vyz*dy*dz)*idr*idr + vrx = vxx*dx + vxy*dy + vxz*dz + vry = vxy*dx + vyy*dy + vyz*dz + vrz = vxz*dx + vyz*dy - (vxx+vyy)*dz + + o[i] += mud * (-2 * vrx + 5 * vrr * dx)*idr5 + o[i+N] += mud * (-2 * vry + 5 * vrr * dy) * idr5 + o[i+2*N] += mud * (-2 * vrz + 5 * vrr * dz) * idr5 + return + + cpdef propulsionR4a( self, double [:] o, double [:] r, double [:] M): + """ + Compute angular velocity due to 4a mode of the slip :math:`v=\pi^{R,4a}\cdot M` + ... + + Parameters + ---------- + o: np.array + An array of angular velocities + An array of size 3*N, + r: np.array + An array of positions + An array of size 3*N, + M: np.array + An array of 4a mode of the slip + An array of size 7*N, + """ + + + cdef int N = self.N, i, j + cdef double dx, dy, dz, idr, idr7, idr9, mrrr, mrrx, mrry, mrrz, mxxx, myyy, mxxy, mxxz, mxyy, mxyz, myyz + cdef double mud = 363 / 80 * self.a ** 4, L = self.L + cdef int neighbors[2] + + for i in prange(N, nogil=True): + neighbors[0] = i - 1 + neighbors[1] = i + 1 + for j in neighbors: + if i == 0 and j == i - 1: + j = N - 1 + dx = r[i] - (r[j] - L) + elif i == N - 1 and j == i + 1: + j = 0 + dx = r[i] - (r[j] + L) + else: + dx = r[i] - r[j] + + mxxx = M[j] + myyy = M[j+N] + mxxy = M[j+2*N] + mxxz = M[j+3*N] + mxyy = M[j+4*N] + mxyz = M[j+5*N] + myyz = M[j+6*N] + + dy = r[i+N] - r[j+N] + dz = r[i+2*N] - r[j+2*N] + idr = 1.0/sqrt( dx*dx + dy*dy + dz*dz ) + idr7 = idr*idr*idr*idr*idr*idr*idr + idr9 = idr7*idr*idr + + mrrr = mxxx*dx*(dx*dx-3*dz*dz) + 3*mxxy*dy*(dx*dx-dz*dz) + mxxz*dz*(3*dx*dx-dz*dz) +\ + 3*mxyy*dx*(dy*dy-dz*dz) + 6*mxyz*dx*dy*dz + myyy*dy*(dy*dy-3*dz*dz) + myyz*dz*(3*dy*dy-dz*dz) + mrrx = mxxx*(dx*dx-dz*dz) + mxyy*(dy*dy-dz*dz) + 2*mxxy*dx*dy + 2*mxxz*dx*dz + 2*mxyz*dy*dz + mrry = mxxy*(dx*dx-dz*dz) + myyy*(dy*dy-dz*dz) + 2*mxyy*dx*dy + 2*mxyz*dx*dz + 2*myyz*dy*dz + mrrz = mxxz*(dx*dx-dz*dz) + myyz*(dy*dy-dz*dz) + 2*mxyz*dx*dy - 2*(mxxx+mxyy)*dx*dz - 2*(mxxy+myyy)*dy*dz + + o[i] += mud * (15 * mrrx * idr7 - 35 * mrrr * dx * idr9) + o[i+N] += mud * (15 * mrry * idr7 - 35 * mrrr * dy * idr9) + o[i+2*N] += mud * (15 * mrrz * idr7 - 35 * mrrr * dz * idr9) + return \ No newline at end of file diff --git a/pystokes/unbounded.pyx b/pystokes/unbounded.pyx index f99a6c2..2721e40 100644 --- a/pystokes/unbounded.pyx +++ b/pystokes/unbounded.pyx @@ -282,6 +282,7 @@ cdef class Rbm: cdef int N = self.N, i, j, Z=2*N cdef double dx, dy, dz, idr, idr5, vxx, vyy, vxy, vxz, vyz, vrx, vry, vrz + cdef double mud = 13.0*self.b*self.b*self.b/12 for i in prange(N, nogil=True): for j in range(N): @@ -300,9 +301,9 @@ cdef class Rbm: vry = vxy*dx + vyy*dy + vyz*dz vrz = vxz*dx + vyz*dy - (vxx+vyy)*dz - v[i] -= 8*( dy*vrz - dz*vry )*idr5 - v[i+N] -= 8*( dz*vrx - dx*vrz )*idr5 - v[i+Z] -= 8*( dx*vry - dy*vrx )*idr5 + v[i] += mud * (vry * dz - vrz * dy) * idr5 + v[i+N] += mud * (vrz * dx - vrx * dz) * idr5 + v[i+2*N] += mud * (vrx * dy - vry * dx) * idr5 return @@ -379,6 +380,7 @@ cdef class Rbm: cdef int N = self.N, i, j cdef double dx, dy, dz, idr, idr7 cdef double mrrx, mrry, mrrz, mxxx, myyy, mxxy, mxxz, mxyy, mxyz, myyz + cdef double mud = -363/8 * self.b ** 4 for i in prange(N, nogil=True): for j in range(N): @@ -399,9 +401,9 @@ cdef class Rbm: mrry = mxxy*(dx*dx-dz*dz) + myyy*(dy*dy-dz*dz) + 2*mxyy*dx*dy + 2*mxyz*dx*dz + 2*myyz*dy*dz mrrz = mxxz*(dx*dx-dz*dz) + myyz*(dy*dy-dz*dz) + 2*mxyz*dx*dy - 2*(mxxx+mxyy)*dx*dz - 2*(mxxy+myyy)*dy*dz - v[i] -= 6*( dy*mrrz - dz*mrry )*idr7 - v[i+N] -= 6*( dz*mrrx - dx*mrrz )*idr7 - v[i+2*N] -= 6*( dx*mrry - dy*mrrx )*idr7 + v[i] -= mud*( dy*mrrz - dz*mrry )*idr7 + v[i+N] -= mud*( dz*mrrx - dx*mrrz )*idr7 + v[i+2*N] -= mud*( dx*mrry - dy*mrrx )*idr7 return @@ -537,7 +539,7 @@ cdef class Rbm: oz += 3*(srx*dy - sry*dx )*idr5 o[i] += ox*mus - o[i+N] += oy*mus + o[i+N] += oy*mus o[i+Z] += oz*mus return @@ -562,6 +564,7 @@ cdef class Rbm: cdef int N = self.N, i, j, Z=2*N cdef double dx, dy, dz, idr, idr2, idr5, vxx, vyy, vxy, vxz, vyz, vrr, vrx, vry, vrz + cdef double mud = 13.0 / 24.0 * self.b * self.b * self.b for i in prange(N, nogil=True): for j in range(N): @@ -581,9 +584,9 @@ cdef class Rbm: vry = vxy*dx + vyy*dy + vyz*dz vrz = vxz*dx + vyz*dy - (vxx+vyy)*dz - o[i] += ( 32*vrx- 20*vrr*dx )*idr5 - o[i+N] += ( 32*vry- 20*vrr*dy )*idr5 - o[i+Z] += ( 32*vrz- 20*vrr*dz )*idr5 + o[i] += mud * (-2 * vrx + 5 * vrr * dx)*idr5 + o[i+N] += mud * (-2 * vry + 5 * vrr * dy) * idr5 + o[i+2*N] += mud * (-2 * vrz + 5 * vrr * dz) * idr5 return @@ -658,6 +661,7 @@ cdef class Rbm: cdef int N = self.N, i, j cdef double dx, dy, dz, idr, idr7, idr9, mrrr, mrrx, mrry, mrrz, mxxx, myyy, mxxy, mxxz, mxyy, mxyz, myyz + cdef double mud = 363 / 80 * self.b ** 4 for i in prange(N, nogil=True): for j in range(N): @@ -682,11 +686,9 @@ cdef class Rbm: mrry = mxxy*(dx*dx-dz*dz) + myyy*(dy*dy-dz*dz) + 2*mxyy*dx*dy + 2*mxyz*dx*dz + 2*myyz*dy*dz mrrz = mxxz*(dx*dx-dz*dz) + myyz*(dy*dy-dz*dz) + 2*mxyz*dx*dy - 2*(mxxx+mxyy)*dx*dz - 2*(mxxy+myyy)*dy*dz - o[i] += 21*mrrr*dx*idr9 - 9*mrrx*idr7 - o[i+N] += 21*mrrr*dy*idr9 - 9*mrry*idr7 - o[i+2*N] += 21*mrrr*dz*idr9 - 9*mrrz*idr7 - else: - pass + o[i] += mud * (15 * mrrx * idr7 - 35 * mrrr * dx * idr9) + o[i+N] += mud * (15 * mrry * idr7 - 35 * mrrr * dy * idr9) + o[i+2*N] += mud * (15 * mrrz * idr7 - 35 * mrrr * dz * idr9) return diff --git a/pystokes/utils.pyx b/pystokes/utils.pyx index 509f2ce..c99fb2a 100644 --- a/pystokes/utils.pyx +++ b/pystokes/utils.pyx @@ -117,12 +117,12 @@ cpdef irreducibleTensors(l, p, Y0=1): if l==3: for i in range(N): - YY[i] = Y0*(p[i]*p[i]*p[i] - 3/5*p[i]); + YY[i] = Y0*(p[i]*p[i]*p[i] - 3/5*p[i]); YY[i+N] = Y0*(p[i+N]*p[i+N]*p[i+N] - 3/5*p[i+N]); YY[i+2*N] = Y0*(p[i]*p[i]*p[i+N] - 1/5*p[i+N]); YY[i+3*N] = Y0*(p[i]*p[i]*p[i+2*N] - 1/5*p[i+2*N]); - YY[i+4*N] = Y0*(p[i]*p[i+N]*p[1+N] -1/5* p[i]); - YY[i+5*N] = Y0*(p[i+N]*p[i+N]*p[i+2*N]); + YY[i+4*N] = Y0*(p[i]*p[i+N]*p[i+N] -1/5* p[i]); + YY[i+5*N] = Y0*(p[i]*p[i+N]*p[i+2*N]); YY[i+6*N] = Y0*(p[i+N]*p[i+N]*p[i+2*N] -1/5*p[i+2*N]); return YY @@ -222,6 +222,72 @@ def simulate(rp0, Tf, Nts, rhs, integrator='odeint', filename='this.mat', savemat(filename, {'X':X, 't':time_points}) return +def simulate_vo2rp(rp0, Tf, Nts, vo, filename='this.mat', + Ti=0, maxNumSteps=100000, **kwargs): + """ + Simulates using choice of integrator + ... + + Parameters + ---------- + rp0 : np.array + Initial condition + Tf : int + Final time + Nts: int + Number of points to return data + vo : Python Function + velcotity and angular velocity + filename: string + filename to write the data. + Deafult is 'this.mat' + """ + + def dxdtEval(rp, t): + """ + returns velocity and angular velocity + """ + return vo(rp) + + def dot(A, b): + """ + A: Array size of (3, 3) + b: Array size of 3 + returns A \dot b + """ + retval = np.zeros(3) + for i in range(3): + for j in range(3): + retval[i] += A[i, j] * b[j] + return retval + + from scipy.integrate import solve_ivp + from scipy.linalg import expm + + time_points=np.linspace(Ti, Tf, Nts+1) + N = int(len(rp0) / 6) + X = np.zeros((Nts + 1, len(rp0))) + X[0] = rp0 + for n in range(Nts): + X[n + 1] = solve_ivp(lambda t, xt: dxdtEval(xt,t), + [time_points[n], time_points[n + 1]], X[n], + t_eval=[time_points[n], time_points[n + 1]], **kwargs).y.T[1] + P = X[n + 1, 3*N:6*N] + P_old = X[n, 3*N:6*N] + O = vo(X[n])[3*N:6*N] + for p in range(N): + Op = O[p::N] + skewOdt = np.array([[0.0, -Op[2], Op[1]], + [Op[2], 0.0, -Op[0]], + [-Op[1], Op[0], 0.0]]) \ + * (time_points[n + 1] - time_points[n]) + P[p::N] = dot(expm(skewOdt), P_old[p::N]) + X[n + 1, 3*N:6*N] = P + + from scipy.io import savemat + savemat(filename, {'X':X, 't':time_points}) + return + @@ -849,4 +915,3 @@ cpdef MSD3d(double [:] d2, double [:] xt, double [:] yt, double [:] zt): dr2 += dx*dx + dy*dy + dz*dz d2[i] = dr2/(Nt-i) return -1