From 5ab0e63bb05a384fd95813d1ad335339e7308f21 Mon Sep 17 00:00:00 2001 From: Takuya Date: Tue, 12 Dec 2023 16:05:19 +0900 Subject: [PATCH 1/6] Update the orientation calculation - Use matrix exponential to calculate the orientation --- pystokes/utils.pyx | 67 ++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 67 insertions(+) diff --git a/pystokes/utils.pyx b/pystokes/utils.pyx index 8f02fea..93e8146 100644 --- a/pystokes/utils.pyx +++ b/pystokes/utils.pyx @@ -221,6 +221,73 @@ 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 + From abf57d15e7839e1e36e99f14426515d40cd0f05d Mon Sep 17 00:00:00 2001 From: Takuya Kobayashi Date: Wed, 22 Jan 2025 11:41:33 +0900 Subject: [PATCH 2/6] implement nearest-neighbor pair approximation --- pystokes/__init__.py | 1 + pystokes/pair.pxd | 35 ++++ pystokes/pair.pyx | 377 +++++++++++++++++++++++++++++++++++++++++++ 3 files changed, 413 insertions(+) create mode 100644 pystokes/pair.pxd create mode 100644 pystokes/pair.pyx diff --git a/pystokes/__init__.py b/pystokes/__init__.py index 7ca0cb6..736240c 100644 --- a/pystokes/__init__.py +++ b/pystokes/__init__.py @@ -1,5 +1,6 @@ import pystokes.interface import pystokes.periodic +import pystokes.pair import pystokes.twoWalls import pystokes.unbounded import pystokes.utils diff --git a/pystokes/pair.pxd b/pystokes/pair.pxd new file mode 100644 index 0000000..cf67507 --- /dev/null +++ b/pystokes/pair.pxd @@ -0,0 +1,35 @@ +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) + + + ## 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) \ No newline at end of file diff --git a/pystokes/pair.pyx b/pystokes/pair.pyx new file mode 100644 index 0000000..6863049 --- /dev/null +++ b/pystokes/pair.pyx @@ -0,0 +1,377 @@ +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 + + + ## 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 \ No newline at end of file From 759ed940d47847d3bda2f8ee9beeef9fe293cdfc Mon Sep 17 00:00:00 2001 From: Takuya Kobayashi Date: Fri, 7 Mar 2025 17:42:37 +0900 Subject: [PATCH 3/6] fixed --- pystokes/unbounded.pyx | 34 +++++++++++++++++++--------------- 1 file changed, 19 insertions(+), 15 deletions(-) diff --git a/pystokes/unbounded.pyx b/pystokes/unbounded.pyx index 7289191..45a3f82 100644 --- a/pystokes/unbounded.pyx +++ b/pystokes/unbounded.pyx @@ -181,7 +181,7 @@ cdef class Rbm: cdef double dx, dy, dz, dr, idr, idr3 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 double sxx, sxy, sxz, syz, syy, srr, srx, sry, srz, mus = (28.0*self.a*self.a)/24 for i in prange(N, nogil=True): vx=0; vy=0; vz=0; @@ -279,6 +279,7 @@ cdef class Rbm: 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 for i in prange(N, nogil=True): for j in range(N): @@ -297,9 +298,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+2*N] -= 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 else: pass 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.a ** 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 else: pass return @@ -510,7 +512,7 @@ cdef class Rbm: cdef int N = self.N, i, j, xx=2*N cdef double dx, dy, dz, idr, idr5, ox, oy, oz - cdef double sxx, sxy, sxz, syz, syy, srr, srx, sry, srz, mus = -(28.0*self.a*self.a)/24 + cdef double sxx, sxy, sxz, syz, syy, srr, srx, sry, srz, mus = (28.0*self.a*self.a)/24 for i in prange(N, nogil=True): ox=0; oy=0; oz=0; @@ -535,7 +537,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+xx] += oz*mus return @@ -560,6 +562,7 @@ cdef class Rbm: 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 for i in prange(N, nogil=True): for j in range(N): @@ -579,9 +582,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+2*N] += ( 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 else : pass 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.a ** 4 for i in prange(N, nogil=True): for j in range(N): @@ -682,9 +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 + 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) else: pass return From 7c7148dea99571abab04d25e48dbe9ddc01dea37 Mon Sep 17 00:00:00 2001 From: Takuya Kobayashi Date: Mon, 14 Apr 2025 13:40:43 +0900 Subject: [PATCH 4/6] implement nearest-neighbor approximations for 1d Cosserat solids - nearest_neighbors.pyx for 1d nearest-neighbor approximations - forceField.pyx for 1d Cosserat solids --- pystokes/__init__.py | 2 +- pystokes/forceFields.pxd | 4 + pystokes/forceFields.pyx | 146 ++++++++++ pystokes/{pair.pxd => nearest_neighbors.pxd} | 16 +- pystokes/{pair.pyx => nearest_neighbors.pyx} | 288 ++++++++++++++++++- 5 files changed, 447 insertions(+), 9 deletions(-) rename pystokes/{pair.pxd => nearest_neighbors.pxd} (68%) rename pystokes/{pair.pyx => nearest_neighbors.pyx} (54%) diff --git a/pystokes/__init__.py b/pystokes/__init__.py index 736240c..421c8ca 100644 --- a/pystokes/__init__.py +++ b/pystokes/__init__.py @@ -1,6 +1,6 @@ import pystokes.interface import pystokes.periodic -import pystokes.pair +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 465b6b8..a4f5ba8 100644 --- a/pystokes/forceFields.pyx +++ b/pystokes/forceFields.pyx @@ -977,7 +977,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) @@ -1060,3 +1115,94 @@ cdef class Torques: T[i+xx] += 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/pair.pxd b/pystokes/nearest_neighbors.pxd similarity index 68% rename from pystokes/pair.pxd rename to pystokes/nearest_neighbors.pxd index cf67507..8dd2a0b 100644 --- a/pystokes/pair.pxd +++ b/pystokes/nearest_neighbors.pxd @@ -16,20 +16,26 @@ cdef class Rbm: 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) \ No newline at end of file + 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/pair.pyx b/pystokes/nearest_neighbors.pyx similarity index 54% rename from pystokes/pair.pyx rename to pystokes/nearest_neighbors.pyx index 6863049..598e44a 100644 --- a/pystokes/pair.pyx +++ b/pystokes/nearest_neighbors.pyx @@ -172,7 +172,7 @@ cdef class Rbm: 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 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): @@ -215,6 +215,170 @@ cdef class Rbm: 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 @@ -338,7 +502,7 @@ cdef class Rbm: 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 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): @@ -374,4 +538,122 @@ cdef class Rbm: o[i] += ox*mus o[i+N] += oy*mus o[i+xx] += oz*mus - return \ No newline at end of file + 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 From 011566b14e800dc2c6320c0c8cef50240bc86125 Mon Sep 17 00:00:00 2001 From: Takuya Kobayashi Date: Mon, 14 Apr 2025 14:25:24 +0900 Subject: [PATCH 5/6] fixed the tab-space inconsistency error --- pystokes/utils.pyx | 126 ++++++++++++++++++++++----------------------- 1 file changed, 62 insertions(+), 64 deletions(-) diff --git a/pystokes/utils.pyx b/pystokes/utils.pyx index b44ccec..a91939a 100644 --- a/pystokes/utils.pyx +++ b/pystokes/utils.pyx @@ -223,71 +223,70 @@ def simulate(rp0, Tf, Nts, rhs, integrator='odeint', filename='this.mat', 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]], + 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 + * (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 @@ -916,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 From a583fb9d3b07e895e4930adcdb60ba65399351e2 Mon Sep 17 00:00:00 2001 From: Takuya Kobayashi Date: Tue, 22 Apr 2025 10:33:26 +0900 Subject: [PATCH 6/6] fixed 3rd irreducible tensors self.b -> self.a --- pystokes/unbounded.pyx | 8 ++++---- pystokes/utils.pyx | 6 +++--- 2 files changed, 7 insertions(+), 7 deletions(-) diff --git a/pystokes/unbounded.pyx b/pystokes/unbounded.pyx index 2d3f104..2721e40 100644 --- a/pystokes/unbounded.pyx +++ b/pystokes/unbounded.pyx @@ -282,7 +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.a*self.a*self.a/12 + cdef double mud = 13.0*self.b*self.b*self.b/12 for i in prange(N, nogil=True): for j in range(N): @@ -380,7 +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.a ** 4 + cdef double mud = -363/8 * self.b ** 4 for i in prange(N, nogil=True): for j in range(N): @@ -564,7 +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.a * self.a * self.a + 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): @@ -661,7 +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.a ** 4 + cdef double mud = 363 / 80 * self.b ** 4 for i in prange(N, nogil=True): for j in range(N): diff --git a/pystokes/utils.pyx b/pystokes/utils.pyx index a91939a..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