Skip to content

Commit

Permalink
self.a -> self.b
Browse files Browse the repository at this point in the history
  • Loading branch information
@rajeshrinet
rajeshrinet committed Sep 28, 2024
1 parent ff3c35d commit 4287b41
Show file tree
Hide file tree
Showing 11 changed files with 156 additions and 150 deletions.
4 changes: 2 additions & 2 deletions examples/ex01-unboundedFlow.ipynb
View file

Large diffs are not rendered by default.

4 changes: 2 additions & 2 deletions pystokes/interface.pxd
View file
Original file line number Diff line number Diff line change
Expand Up @@ -13,7 +13,7 @@ cdef class Rbm:

cdef readonly int Nx, Ny, Nz, N
cdef readonly np.ndarray fkx, fky, fkz, vkx, vky, vkz, fk0, fx0, Mobility
cdef readonly double Lx, Ly, Lz, a, facx, facy, facz, eta, mu, muv, mur
cdef readonly double Lx, Ly, Lz, b, facx, facy, facz, eta, mu, muv, mur


cpdef mobilityTT(self, double [:] v, double [:] r, double [:] F, double ll=?)
Expand Down Expand Up @@ -61,7 +61,7 @@ cdef class Rbm:
@cython.cdivision(True)
@cython.nonecheck(False)
cdef class Flow:
cdef double a, eta
cdef double b, eta
cdef int N
cdef int Nt

Expand Down
65 changes: 33 additions & 32 deletions pystokes/interface.pyx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -31,15 +31,16 @@ cdef class Rbm:


def __init__(self, radius=1, particles=1, viscosity=1.0):
self.a = radius
self.b = radius
self.N = particles
self.eta = viscosity
self.mu = 1.0/(6*PI*self.eta*self.a)
self.mu = 1.0/(6*PI*self.eta*self.b)
self.muv = 1.0/(8*PI*self.eta)
self.mur = 1.0/(8*PI*self.eta*self.a**3)
self.mur = 1.0/(8*PI*self.eta*self.b**3)

self.Mobility = np.zeros( (3*self.N, 3*self.N), dtype=np.float64)


cpdef mobilityTT(self, double [:] v, double [:] r, double [:] F, double ll=0):
"""
Compute velocity due to body forces using :math:`v=\mu^{TT}\cdot F`
Expand Down Expand Up @@ -92,11 +93,11 @@ cdef class Rbm:
cdef int i, j, N=self.N, xx=2*N
cdef double dx, dy, dz, idr, idr3, idr5, Fdotidr, h2, hsq, tempF
cdef double vx, vy, vz, F1, F2, F3
cdef double mu=self.mu, muv=self.muv, a2=self.a*self.a/3.0
cdef double mu=self.mu, muv=self.muv, a2=self.b*self.b/3.0
cdef double ll1 = (1-ll)/(1+ll), ll2 = ll/(1+ll);

cdef double llp = 1.0/(1+ll);
cdef double a = self.a
cdef double a = self.b
cdef double h, hbar_inv, hbar_inv3, hbar_inv5
cdef double muTTpara1 = 3*(2-3*ll)*llp/16.0, muTTpara2 = (1+2*ll)*llp/16.0
cdef double muTTpara3 = - ll*llp/16.0
Expand All @@ -117,7 +118,7 @@ cdef class Rbm:
Fdotidr = (F[j] * dx + F[j+N] * dy + F[j+xx] * dz)*idr*idr
#
vx += (F[j] +Fdotidr*dx)*idr + a2*(2*F[j] -6*Fdotidr*dx)*idr3
vy += (F[j+N]+Fdotidr*dy)*idr + a2*(2*F[j+N]-6*Fdotidr*dy)*idr3
vy += (F[j+N]+Fdotidr*dy)*idr + a2*(2*F[j+N]-6*Fdotidr*dy)*idr3
vz += (F[j+xx]+Fdotidr*dz)*idr + a2*(2*F[j+xx]-6*Fdotidr*dz)*idr3

##contributions from the image
Expand Down Expand Up @@ -195,9 +196,9 @@ cdef class Rbm:
cdef double ll1 = (1-ll)/(1+ll), ll2 = ll/(1+ll);

cdef double llp = 1.0/(1+ll);
cdef double a = self.a
cdef double a = self.b
cdef double h, hbar_inv, hbar_inv2, hbar_inv4
cdef double muTR0 = 4.0/(3*self.a*self.a)
cdef double muTR0 = 4.0/(3*self.b*self.b)
cdef double muTR1 = -3*llp/16.0, muTR2 = 3*ll2/32.0
cdef double muTR

Expand Down Expand Up @@ -275,11 +276,11 @@ cdef class Rbm:
cdef double dx, dy, dz, idr, idr2, idr3, idr5, idr7, aidr2, trS, h2, hsq
cdef double sxx, syy, szz, sxy, syx, syz, szy, sxz, szx, srr, srx, sry, srz
cdef double Sljrlx, Sljrly, Sljrlz, Sljrjx, Sljrjy, Sljrjz
cdef double vx, vy, vz, mus =(28.0*self.a**3)/24
cdef double vx, vy, vz, mus =(28.0*self.b**3)/24
cdef double ll2 = ll/(1+ll);

cdef double llp = 1.0/(1+ll);
cdef double a = self.a
cdef double a = self.b
cdef double h, hbar_inv, hbar_inv2, hbar_inv4, hbar_inv6
cdef double piT2s11 = 5*ll2/16.0, piT2s12 = -(1+3*ll)*llp/12.0
cdef double piT2s13 = 5*ll2/48.0
Expand Down Expand Up @@ -399,11 +400,11 @@ cdef class Rbm:

cdef int N=self.N, i, j, xx=2*N
cdef double dx, dy, dz, idr, idr3, idr5, Ddotidr, tempD, hsq, h2, D1, D2, D3
cdef double vx, vy, vz, mud = 3.0*self.a*self.a*self.a/5, muv = -1.0*(self.a**5)/10
cdef double vx, vy, vz, mud = 3.0*self.b*self.b*self.b/5, muv = -1.0*(self.b**5)/10
cdef double ll1 = (1-ll)/(1+ll), ll2 = ll/(1+ll);

cdef double llp = 1.0/(1+ll);
cdef double a = self.a
cdef double a = self.b
cdef double h, hbar_inv, hbar_inv3, hbar_inv5
cdef double piT3tpara1 = -(1+2*ll)*llp/80.0
cdef double piT3tpara2 = ll2/40.0
Expand Down Expand Up @@ -478,9 +479,9 @@ cdef class Rbm:
cdef double ll1 = (1-ll)/(1+ll), ll2 = ll/(1+ll)

cdef double llp = 1.0/(1+ll);
cdef double a = self.a
cdef double a = self.b
cdef double h, hbar_inv, hbar_inv2, hbar_inv4
cdef double muRT0 = 4.0/(3*self.a*self.a)
cdef double muRT0 = 4.0/(3*self.b*self.b)
cdef double muRT1 = 3*llp/16.0, muRT2 = -3*ll2/32.0
cdef double muRT, F1, F2

Expand Down Expand Up @@ -540,7 +541,7 @@ cdef class Rbm:
cdef double ox, oy, oz, mur=self.mur, muv=self.muv

cdef double llp = 1.0/(1+ll);
cdef double a = self.a
cdef double a = self.b
cdef double h, hbar_inv, hbar_inv3
cdef double muRRpara = (1-5*ll)*llp/16.0
cdef double muRRperp = (1-ll)*llp/8.0
Expand Down Expand Up @@ -623,7 +624,7 @@ cdef class Rbm:
cdef double ll1 = (1-ll)/(1+ll), ll2 = ll/(1+ll);

cdef double llp = 1.0/(1+ll);
cdef double a = self.a
cdef double a = self.b
cdef double h, hbar_inv, hbar_inv4
cdef double piR3t0 = 3*ll2/(80.0*a)
cdef double piR3t, V1, V2
Expand Down Expand Up @@ -688,7 +689,7 @@ cdef class Rbm:
cdef double ll2 = ll/(1+ll);

cdef double llp = 1.0/(1+ll);
cdef double a = self.a, a_inv = 1.0/a
cdef double a = self.b, a_inv = 1.0/a
cdef double h, hbar_inv, hbar_inv3, hbar_inv5
cdef double piR2s1 = 5.0/32.0, piR2s2 = -ll2/8.0
cdef double piR2s
Expand Down Expand Up @@ -746,8 +747,8 @@ cdef class Rbm:
"""
cdef int i, j, N=self.N, xx=2*N
cdef double dx, dy, dz, idr, h2, hsq, idr2, idr3, idr4, idr5
cdef double mu=self.mu, muv=2*mu*self.a*0.75, a2=self.a*self.a/3.0
cdef double vx, vy, vz, mm=1.0/(.75*self.a)
cdef double mu=self.mu, muv=2*mu*self.b*0.75, a2=self.b*self.b/3.0
cdef double vx, vy, vz, mm=1.0/(.75*self.b)

cdef double [:, :] M = self.Mobility
cdef double [:] Fr = np.random.normal(size=3*N)
Expand Down Expand Up @@ -851,7 +852,7 @@ cdef class Rbm:

cdef int i, j, N=self.N, xx=2*N
cdef double dx, dy, dz, idr, h2, hsq, idr2, idr3, idr4, idr5
cdef double mur=self.muv, muv=0.25*sqrt(2.0)*mur, mm=4.0/(self.a**3)
cdef double mur=self.muv, muv=0.25*sqrt(2.0)*mur, mm=4.0/(self.b**3)
cdef double ox, oy, oz

cdef double [:, :] M = self.Mobility
Expand Down Expand Up @@ -954,14 +955,14 @@ cdef class Rbm:

cdef double [:] Fr = np.random.normal(size=3*N)

cdef double a = self.a
cdef double a = self.b
cdef double h, hbar_inv, hbar_inv2, hbar_inv3, hbar_inv4, hbar_inv5
cdef double muTTparaCoeff1 = 3*(2-3*ll)*llp/16.0, muTTparaCoeff2 = (1+2*ll)*llp/16.0
cdef double muTTparaCoeff3 = - ll*llp/16.0
cdef double muTTperpCoeff1 = - 3*(2+3*ll)*llp/8.0, muTTperpCoeff2 = (1+4*ll)*llp/8.0
cdef double muTTperpCoeff3 = - ll*llp/8.0

cdef double muTRCoeff = 4.0/(3*self.a*self.a)
cdef double muTRCoeff = 4.0/(3*self.b*self.b)
cdef double muTRCoeff1 = -3*llp/16.0, muTRCoeff2 = 3*ll2/32.0

cdef double muRRparaCoeff = (1-5*ll)*llp/16.0
Expand Down Expand Up @@ -1023,14 +1024,14 @@ cdef class Rbm:

cdef double [:] Tr = np.random.normal(size=3*N)

cdef double a = self.a
cdef double a = self.b
cdef double h, hbar_inv, hbar_inv2, hbar_inv3, hbar_inv4, hbar_inv5
cdef double muTTparaCoeff1 = 3*(2-3*ll)*llp/16.0, muTTparaCoeff2 = (1+2*ll)*llp/16.0
cdef double muTTparaCoeff3 = - ll*llp/16.0
cdef double muTTperpCoeff1 = - 3*(2+3*ll)*llp/8.0, muTTperpCoeff2 = (1+4*ll)*llp/8.0
cdef double muTTperpCoeff3 = - ll*llp/8.0

cdef double muTRCoeff = 4.0/(3*self.a*self.a)
cdef double muTRCoeff = 4.0/(3*self.b*self.b)
cdef double muTRCoeff1 = -3*llp/16.0, muTRCoeff2 = 3*ll2/32.0

cdef double muRRparaCoeff = (1-5*ll)*llp/16.0
Expand Down Expand Up @@ -1088,14 +1089,14 @@ cdef class Rbm:

cdef double [:] Fr = np.random.normal(size=3*N)

cdef double a = self.a
cdef double a = self.b
cdef double h, hbar_inv, hbar_inv2, hbar_inv3, hbar_inv4, hbar_inv5
cdef double muTTparaCoeff1 = 3*(2-3*ll)*llp/16.0, muTTparaCoeff2 = (1+2*ll)*llp/16.0
cdef double muTTparaCoeff3 = - ll*llp/16.0
cdef double muTTperpCoeff1 = - 3*(2+3*ll)*llp/8.0, muTTperpCoeff2 = (1+4*ll)*llp/8.0
cdef double muTTperpCoeff3 = - ll*llp/8.0

cdef double muTRCoeff = 4.0/(3*self.a*self.a)
cdef double muTRCoeff = 4.0/(3*self.b*self.b)
cdef double muTRCoeff1 = -3*llp/16.0, muTRCoeff2 = 3*ll2/32.0

cdef double muRRparaCoeff = (1-5*ll)*llp/16.0
Expand Down Expand Up @@ -1153,14 +1154,14 @@ cdef class Rbm:

cdef double [:] Tr = np.random.normal(size=3*N)

cdef double a = self.a
cdef double a = self.b
cdef double h, hbar_inv, hbar_inv2, hbar_inv3, hbar_inv4, hbar_inv5
cdef double muTTparaCoeff1 = 3*(2-3*ll)*llp/16.0, muTTparaCoeff2 = (1+2*ll)*llp/16.0
cdef double muTTparaCoeff3 = - ll*llp/16.0
cdef double muTTperpCoeff1 = - 3*(2+3*ll)*llp/8.0, muTTperpCoeff2 = (1+4*ll)*llp/8.0
cdef double muTTperpCoeff3 = - ll*llp/8.0

cdef double muTRCoeff = 4.0/(3*self.a*self.a)
cdef double muTRCoeff = 4.0/(3*self.b*self.b)
cdef double muTRCoeff1 = -3*llp/16.0, muTRCoeff2 = 3*ll2/32.0

cdef double muRRparaCoeff = (1-5*ll)*llp/16.0
Expand Down Expand Up @@ -1237,7 +1238,7 @@ cdef class Flow:
"""

def __init__(self, radius=1, particles=1, viscosity=1, gridpoints=32):
self.a = radius
self.b = radius
self.N = particles
self.Nt = gridpoints
self.eta= viscosity
Expand Down Expand Up @@ -1290,7 +1291,7 @@ cdef class Flow:

cdef int i, j, N=self.N, Nt=self.Nt, xx=2*N
cdef double dx, dy, dz, idr, idr3, idr5, Fdotidr, tempF, hsq, h2, F3
cdef double vx, vy, vz, muv=1.0/(8*PI*self.eta), a2=self.a*self.a/6.0
cdef double vx, vy, vz, muv=1.0/(8*PI*self.eta), a2=self.b*self.b/6.0

for i in prange(Nt, nogil=True):
vx=0; vy=0; vz=0;
Expand Down Expand Up @@ -1472,7 +1473,7 @@ cdef class Flow:
cdef double dx, dy, dz, idr, idr2, idr3, idr5, idr7, aidr2, trS, h2, hsq
cdef double sxx, syy, szz, sxy, syx, syz, szy, sxz, szx, srr, srx, sry, srz
cdef double Sljrlx, Sljrly, Sljrlz, Sljrjx, Sljrjy, Sljrjz
cdef double vx, vy, vz, mus = (28.0*self.a**3)/24
cdef double vx, vy, vz, mus = (28.0*self.b**3)/24

for i in prange(Nt, nogil=True):
vx=0; vy=0; vz=0;
Expand Down Expand Up @@ -1571,7 +1572,7 @@ cdef class Flow:
"""
cdef int i, j, N=self.N, Nt=self.Nt, xx=2*N
cdef double dx, dy, dz, idr, idr3, idr5, Ddotidr, tempD, hsq, h2, D3
cdef double vx, vy, vz, mud = 3.0*self.a*self.a*self.a/5, muv = -1.0*(self.a**5)/10
cdef double vx, vy, vz, mud = 3.0*self.b*self.b*self.b/5, muv = -1.0*(self.b**5)/10

for i in prange(Nt, nogil=True):
vx=0; vy=0; vz=0;
Expand Down
6 changes: 3 additions & 3 deletions pystokes/periodic.pxd
View file
Original file line number Diff line number Diff line change
Expand Up @@ -8,7 +8,7 @@ from cython.parallel import prange
@cython.nonecheck(False)
cdef class Rbm:
cdef int N
cdef double a, eta, L, mu, xi, muv, mur
cdef double b, eta, L, mu, xi, muv, mur

cpdef mobilityTT(self, double [:] v, double [:] r, double [:] F, int Nb=?, int Nm=?, double xi0=?)

Expand Down Expand Up @@ -57,9 +57,9 @@ cdef class Rbm:
@cython.cdivision(True)
@cython.nonecheck(False)
cdef class Flow:
cdef double a
cdef double b
cdef int N, Nt
cdef double L, eta
cdef double L, eta, xi

cpdef flowField1s(self, double [:] v, double [:] rt, double [:] r, double [:] F, int Nb=?, int Nm=?, double xi0=?)

Expand Down
Loading

0 comments on commit 4287b41

Please sign in to comment.