Skip to content

Commit

Permalink
rename 'Np' as 'N'
Browse files Browse the repository at this point in the history
  • Loading branch information
@rajeshrinet
rajeshrinet committed Sep 12, 2023
1 parent 945aa61 commit 2265279
Show file tree
Hide file tree
Showing 27 changed files with 1,664 additions and 1,766 deletions.
6 changes: 3 additions & 3 deletions examples/ex1-unboundedFlow.ipynb
View file

Large diffs are not rendered by default.

10 changes: 5 additions & 5 deletions examples/ex2-flowPlaneSurface.ipynb
View file

Large diffs are not rendered by default.

39 changes: 20 additions & 19 deletions examples/ex3-crystalNucleation.ipynb
View file

Large diffs are not rendered by default.

48 changes: 22 additions & 26 deletions examples/ex4-crystallization.ipynb
View file

Large diffs are not rendered by default.

22 changes: 10 additions & 12 deletions examples/ex5-phoreticField.ipynb
View file

Large diffs are not rendered by default.

32 changes: 15 additions & 17 deletions examples/ex6-arrestedCluster.ipynb
View file

Large diffs are not rendered by default.

14 changes: 7 additions & 7 deletions examples/ex7-benchmark.ipynb
View file

Large diffs are not rendered by default.

100 changes: 50 additions & 50 deletions examples/ex8-holographicTrap.ipynb
View file
Original file line number Diff line number Diff line change
Expand Up @@ -44,9 +44,9 @@
"source": [
"DTYPE = np.float64\n",
"class trap:\n",
" def __init__(self, a, Np, vs, eta, dim, S0, D0, k, ljeps, ljrmin):\n",
" def __init__(self, a, N, vs, eta, dim, S0, D0, k, ljeps, ljrmin):\n",
" self.a = a\n",
" self.Np = Np\n",
" self.N = N\n",
" self.vs = vs\n",
" self.eta = eta\n",
" self.dim = dim\n",
Expand All @@ -56,9 +56,9 @@
" self.k = k\n",
" self.ljrmin = ljrmin\n",
" self.ljeps = ljeps\n",
" self.dX = np.zeros( 6*self.Np, dtype=DTYPE)\n",
" self.dX = np.zeros( 6*self.N, dtype=DTYPE)\n",
" \n",
" self.uRbm = pystokes.unbounded.Rbm(self.a, self.Np, self.eta)\n",
" self.uRbm = pystokes.unbounded.Rbm(self.a, self.N, self.eta)\n",
"\n",
"\n",
" def initialise(self, initialConfig, trapCentre):\n",
Expand All @@ -67,32 +67,32 @@
" return\n",
"\n",
" def rhs(self, rp):\n",
" Np = self.Np \n",
" xx=2*Np; vs = self.vs\n",
" N = self.N \n",
" xx=2*N; vs = self.vs\n",
" \n",
" r = rp[0:3*Np]\n",
" p = rp[3*Np:6*Np] \n",
" r = rp[0:3*N]\n",
" p = rp[3*N:6*N] \n",
" F = - self.k*(r-trapCentre)\n",
" \n",
" v = r*0; o=r*0;\n",
" self.uRbm.mobilityTT(v, r, F)\n",
" self.uRbm.mobilityRT(o, r, F)\n",
" \n",
" self.dX[0:3*Np] = vs*p + v \n",
" self.dX[3*Np:4*Np] = o[Np:2*Np ]*p[2*Np:3*Np] - o[2*Np:3*Np]*p[Np:2*Np]\n",
" self.dX[4*Np:5*Np] = o[2*Np:3*Np]*p[0:Np ] - o[0:Np ]*p[2*Np:3*Np]\n",
" self.dX[5*Np:6*Np] = o[0:Np ]*p[Np:2*Np ] - o[Np:2*Np ]*p[0:Np ]\n",
" self.dX[0:3*N] = vs*p + v \n",
" self.dX[3*N:4*N] = o[N:2*N ]*p[2*N:3*N] - o[2*N:3*N]*p[N:2*N]\n",
" self.dX[4*N:5*N] = o[2*N:3*N]*p[0:N ] - o[0:N ]*p[2*N:3*N]\n",
" self.dX[5*N:6*N] = o[0:N ]*p[N:2*N ] - o[N:2*N ]*p[0:N ]\n",
" return self.dX\n",
" \n",
" def simulate(self, dt, N):\n",
" '''run simulation and save data'''\n",
" X = np.zeros( (N+1, 6*self.Np), dtype=DTYPE)\n",
" X = np.zeros( (N+1, 6*self.N), dtype=DTYPE)\n",
" X[0, :] = self.rp0\n",
"\n",
" for i in range(N):\n",
" X[i+1, :] = X[i, :] + self.rhs(X[i, :])*dt\n",
" \n",
" sio.savemat('Np=%s_vs=%4.4f_K=%4.4f.mat'%(self.Np, self.vs, self.k), {'trapCentre':self.trapCentre, 'X':X, 't':dt, 'Np':self.Np,'k':self.k, 'vs':self.vs, 'S0':self.S0,})\n",
" sio.savemat('N=%s_vs=%4.4f_K=%4.4f.mat'%(self.N, self.vs, self.k), {'trapCentre':self.trapCentre, 'X':X, 't':dt, 'N':self.N,'k':self.k, 'vs':self.vs, 'S0':self.S0,})\n",
" return"
]
},
Expand All @@ -103,7 +103,7 @@
"outputs": [],
"source": [
"## parameters and initial conditions\n",
"Np = 9 # number of particles\n",
"N = 9 # number of particles\n",
"vs = 1 # self-propulsion speed\n",
"A = 2 # number vs/k\n",
"k = vs/A # stiffness of the trap\n",
Expand All @@ -115,54 +115,54 @@
"\n",
"S0, D0 = 0.01, 0.01 # strength of the stresslet and potDiple\n",
"ljrmin, ljeps = 3, .01 # lennard-jones parameters\n",
"dt, Npts = 10, 200 # final time and number of points \n",
"dt, Nts = 10, 200 # final time and number of points \n",
"\n",
"# instantiate the class trap from trap.pyx for simulation \n",
"rm = trap(b, Np, vs, eta, dim, S0, D0, k, ljeps, ljrmin) \n",
"rm = trap(b, N, vs, eta, dim, S0, D0, k, ljeps, ljrmin) \n",
"\n",
"# set initial condition\n",
"def initialConfig(rp0, trapCentre, theta, b, a0, vs, k, Np): \n",
"def initialConfig(rp0, trapCentre, theta, b, a0, vs, k, N): \n",
" '''method for preparing an initial system''' \n",
" rr = (6*np.pi*eta*vs*b)/k # confinement radius\n",
" t1 = np.pi/180 \n",
" \n",
" if Np==1: \n",
" if N==1: \n",
" rp0[0], rp0[1], rp0[2] = 0, 0, 8 # Position \n",
" rp0[3], rp0[4], rp0[5] = 0, 0, -1 # Orientation \n",
" \n",
" elif Np==3: \n",
" elif N==3: \n",
" trapCentre[0], trapCentre[3] = 0, 0 \n",
" trapCentre[1], trapCentre[4] = a0, a0\n",
" trapCentre[2], trapCentre[5] = -a0, a0 \n",
" theta[0:3] = 90*np.pi/180 \n",
" for i in range(Np):\n",
" for i in range(N):\n",
" rp0[i ] = trapCentre[i ] + rr*np.cos(theta[i])\n",
" rp0[i+Np ] = trapCentre[i+Np] + rr*np.sin(theta[i])\n",
" rp0[i+3*Np] = np.cos(theta[i])\n",
" rp0[i+4*Np] = np.sin(theta[i])\n",
" rp0[i+N ] = trapCentre[i+N] + rr*np.sin(theta[i])\n",
" rp0[i+3*N] = np.cos(theta[i])\n",
" rp0[i+4*N] = np.sin(theta[i])\n",
" else:\n",
" Np1d = np.int32(np.round( (Np)**(1.0/2)))\n",
" nnd = Np1d/2 - 0.5; h0=0\n",
" Np= np.int32(Np1d*Np1d)\n",
" for i in range(Np1d):\n",
" for j in range(Np1d):\n",
" ii = i*Np1d + j\n",
" N1d = np.int32(np.round( (N)**(1.0/2)))\n",
" nnd = N1d/2 - 0.5; h0=0\n",
" N= np.int32(N1d*N1d)\n",
" for i in range(N1d):\n",
" for j in range(N1d):\n",
" ii = i*N1d + j\n",
" trapCentre[ii] = a0*(-nnd + i)\n",
" trapCentre[ii+Np] = a0*(-nnd + j)\n",
" trapCentre[ii+2*Np] = h0\n",
" theta = np.ones(Np)*np.pi/2\n",
" #theta = np.random.random(Np)*np.pi/2\n",
" for i in range(Np):\n",
" trapCentre[ii+N] = a0*(-nnd + j)\n",
" trapCentre[ii+2*N] = h0\n",
" theta = np.ones(N)*np.pi/2\n",
" #theta = np.random.random(N)*np.pi/2\n",
" for i in range(N):\n",
" rp0[i ] = trapCentre[i ] + rr*np.cos(theta[i])\n",
" rp0[i+Np ] = trapCentre[i+Np] + rr*np.sin(theta[i])\n",
" rp0[i+3*Np] = np.cos(theta[i])\n",
" rp0[i+4*Np] = np.sin(theta[i])\n",
" rp0[i+N ] = trapCentre[i+N] + rr*np.sin(theta[i])\n",
" rp0[i+3*N] = np.cos(theta[i])\n",
" rp0[i+4*N] = np.sin(theta[i])\n",
" \n",
"rp0 = np.zeros(6*Np) # memory allocation for positions and orientation of colloids\n",
"trapCentre = np.zeros(3*Np) # memory allocation for trap centers\n",
"theta = np.zeros(Np) # angle of the colloids about the trap centers\n",
"rp0 = np.zeros(6*N) # memory allocation for positions and orientation of colloids\n",
"trapCentre = np.zeros(3*N) # memory allocation for trap centers\n",
"theta = np.zeros(N) # angle of the colloids about the trap centers\n",
"\n",
"initialConfig(rp0, trapCentre, theta, b, a0, vs, k, Np)\n",
"initialConfig(rp0, trapCentre, theta, b, a0, vs, k, N)\n",
"rm.initialise(rp0, trapCentre)"
]
},
Expand All @@ -173,7 +173,7 @@
"outputs": [],
"source": [
"# simulate the system and save data\n",
"rm.simulate(dt, Npts)"
"rm.simulate(dt, Nts)"
]
},
{
Expand All @@ -191,12 +191,12 @@
"source": [
"%matplotlib inline\n",
"\n",
"data = sio.loadmat('Np=9_vs=1.0000_K=0.5000.mat')\n",
"data = sio.loadmat('N=9_vs=1.0000_K=0.5000.mat')\n",
"X = data['X']\n",
"tm = data['t']\n",
"k = data['k']\n",
"vs = data['vs']\n",
"Np = data['Np']\n",
"N = data['N']\n",
"tC = data['trapCentre']"
]
},
Expand All @@ -222,10 +222,10 @@
"\n",
"def plotConfig(n, n_):\n",
" ax = fig.add_subplot(1, 5, n_, aspect='equal', )\n",
" for i in range(int(Np)):\n",
" ax.add_patch(patches.Circle((tC[0, i], tC[0, i+Np]), rr, color='#348abd', alpha=0.32))\n",
" x, y = X[n,i], X[n,Np+i]\n",
" px, py = X[n,3*Np+i], X[n,4*Np+i]\n",
" for i in range(int(N)):\n",
" ax.add_patch(patches.Circle((tC[0, i], tC[0, i+N]), rr, color='#348abd', alpha=0.32))\n",
" x, y = X[n,i], X[n,N+i]\n",
" px, py = X[n,3*N+i], X[n,4*N+i]\n",
" plt.quiver(x,y,px,py)\n",
" plt.title('Time=%d'%n, fontsize=20); plt.axis('off')\n",
"\n",
Expand All @@ -252,7 +252,7 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.9.12"
"version": "3.10.9"
}
},
"nbformat": 4,
Expand Down
2 changes: 1 addition & 1 deletion examples/others/abp/active_brownian_2D.ipynb
View file
Original file line number Diff line number Diff line change
Expand Up @@ -197,7 +197,7 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.9.12"
"version": "3.10.9"
}
},
"nbformat": 4,
Expand Down
30 changes: 12 additions & 18 deletions examples/others/fluid_flow/ex-unboundedFlow.ipynb
View file

Large diffs are not rendered by default.

56 changes: 25 additions & 31 deletions examples/others/forceFields/ex1.ipynb
View file

Large diffs are not rendered by default.

12 changes: 6 additions & 6 deletions examples/others/periodic/compMobility.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,14 +5,14 @@

#Parameters
a, eta, dim = 1.0, 1.0/6, 3
Np, Nb, Nm = 1, 1, 8
N, Nb, Nm = 1, 1, 8
ta =(4*np.pi/3)**(1.0/3)
L = ta/np.arange(0.01, 0.4, 0.01)

#Memory allocation
v = np.zeros(dim*Np)
r = np.zeros(dim*Np)
F = np.zeros(dim*Np)
v = np.zeros(dim*N)
r = np.zeros(dim*N)
F = np.zeros(dim*N)
vv = np.zeros(np.size(L))
phi = np.zeros(np.size(L) )

Expand All @@ -24,10 +24,10 @@

r[0], r[1], r[2] = 0.0, 0.0, 0.0

ff = pystokes.forceFields.Forces(Np)
ff = pystokes.forceFields.Forces(N)
ff.sedimentation(F, g=-1)

pRbm = pystokes.periodic.Rbm(a, Np, eta, L[i])
pRbm = pystokes.periodic.Rbm(a, N, eta, L[i])
pRbm.mobilityTT(v, r, F, Nb, Nm)

phi[i] = (4*np.pi*a**3)/(3*L[i]**3)
Expand Down
64 changes: 32 additions & 32 deletions examples/others/periodic/crowleyInstability.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6,59 +6,59 @@



a, Np1d, Np = 1, 10, 100 # radius and number of particles
a, N1d, N = 1, 10, 100 # radius and number of particles
L, dim = 128, 3 # size and dimensionality and the box
latticeShape = 'square' # shape of the lattice
v = np.zeros(dim*Np) # Memory allocation for velocity
r = np.zeros(dim*Np) # Position vector of the particles
F = np.zeros(dim*Np) # Forces on the particles
v = np.zeros(dim*N) # Memory allocation for velocity
r = np.zeros(dim*N) # Position vector of the particles
F = np.zeros(dim*N) # Forces on the particles
Nb, Nm = 1, 4

rm = pystokes.periodic.Rbm(a, Np, 1.0/6, L) # instantiate the classes
ff = pystokes.forceFields.Forces(Np)
rm = pystokes.periodic.Rbm(a, N, 1.0/6, L) # instantiate the classes
ff = pystokes.forceFields.Forces(N)

def initialise(r, Np1d, shape):
def initialise(r, N1d, shape):
''' this is the module to initialise the particles on the lattice'''
if shape=='cube':
for i in range(Np1d):
for j in range(Np1d):
for k in range(Np1d):
ii = i*Np1d**2 + j*Np1d + k
r[ii] = L/2 + 3*(-Np1d + 2*i)
r[ii+Np] = L/2 + 3*(-Np1d + 2+ 2*j)
r[ii+2*Np] = L/2 + 3*(-Np1d + 2+ 2*k)
for i in range(N1d):
for j in range(N1d):
for k in range(N1d):
ii = i*N1d**2 + j*N1d + k
r[ii] = L/2 + 3*(-N1d + 2*i)
r[ii+N] = L/2 + 3*(-N1d + 2+ 2*j)
r[ii+2*N] = L/2 + 3*(-N1d + 2+ 2*k)
elif shape=='square':
for i in range(Np1d):
for j in range(Np1d):
ii = i*Np1d + j
r[ii] = L/2 + 4*(-Np1d + 2*i)
r[ii+Np] = L/2 + 4*(-Np1d + 2+ 2*j)
for i in range(Np):
r[i+2*Np] = L/2
for i in range(N1d):
for j in range(N1d):
ii = i*N1d + j
r[ii] = L/2 + 4*(-N1d + 2*i)
r[ii+N] = L/2 + 4*(-N1d + 2+ 2*j)
for i in range(N):
r[i+2*N] = L/2

elif shape=='rectangle':
for i in range(Np1d):
for j in range(Np1d):
ii = i*Np1d + j
r[ii] = L/2 + 4*(-Np1d + 2*i)
r[ii+Np] = L/2 + 4*(-Np1d + 2+ 2*j)
for i in range(Np):
r[i+2*Np] = L/2
for i in range(N1d):
for j in range(N1d):
ii = i*N1d + j
r[ii] = L/2 + 4*(-N1d + 2*i)
r[ii+N] = L/2 + 4*(-N1d + 2+ 2*j)
for i in range(N):
r[i+2*N] = L/2
else:
pass


initialise(r, Np1d, latticeShape)
initialise(r, N1d, latticeShape)
dt = 0.01
fig = plt.figure()
for tt in range(32):
ff.sedimentation(F, g=-10) # call the Sedimentation module of ForceFields
v=v*0 # setting v=0 in each time step
rm.mobilityTT(v, r, F, Nb, Nm) # and StokesletV module of pystokes
r = (r + v*dt)%L
x = r[0:Np]
y = r[Np:2*Np]
z = r[2*Np:3*Np]
x = r[0:N]
y = r[N:2*N]
z = r[2*N:3*N]
cc = x*x + y*y + z*z
ax3D = fig.add_subplot(111, projection='3d')
scatCollection = ax3D.scatter(x, y, z, s=30, c=cc, cmap=plt.cm.RdBu )
Expand Down
234 changes: 60 additions & 174 deletions examples/others/periodic/mobilitySedimentingLattice.ipynb
View file

Large diffs are not rendered by default.

5 changes: 3 additions & 2 deletions pystokes/forceFields.pxd
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,7 +7,7 @@ from cython.parallel import prange
@cython.cdivision(True)
@cython.nonecheck(False)
cdef class Forces:
cdef int Np
cdef int N

cpdef VdW(self, double [:] F, double [:] r, double A=?, double a0=?)

Expand Down Expand Up @@ -68,7 +68,8 @@ cdef class Forces:
@cython.cdivision(True)
@cython.nonecheck(False)
cdef class Torques:
cdef int Np
cdef int N

cpdef bottomHeaviness(self, double [:] T, double [:] p, double bh=?)

cpdef magnetic(self, double[:] T, double [:] p, double m0, double Bx, double By, double Bz)

0 comments on commit 2265279

Please sign in to comment.