Simulating magnetic particles in a field free point generated by two opposing magnets

This is probably a long shot with such a short time, but I've been trying to get theoretical data for a project I'm working on. The project involves using a very simplified version of magnetic particle imaging to analyze fluids.

So I've written some code with some help, and I'm pretty sure all the physics is correct but when I run it with time steps (deltaT) of 5E-5 the code calculates velocities of the order of thousands when it should be in millimeters per second. But when the code has time steps of 5E-6 it runs correctly. The problem is the time it takes to complete these simulations at this time step is huge and in my actual experiments, some of these took over 20 minutes.

Basically, I'm wondering if it's possible to either significantly decrease the run time at this time step or be able to increase the time step of the simulation.

Here's my code:

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation

#Arrays for storing time
#t1 = [0]
#t2 = [0]
#t3 = [0]
#t4 = [0]
#t5 = [0]
t6 = [0]

#Time Increment
deltaT = 5E-7
#Number of Particles
N = 100

#Viscosity
mu_w = 0.0010049
mu = 3*mu_w

#Geometry
d = 0.023
L = 0.02

#Constants
mu_0 = 1.25663706E-6

#Particle Properties
m = 1.663223552E-13
moment = 60*m

x, y = np.mgrid[-1.15:1.15:151j, -1:1:151j]*1E-2

#Method for Magnetic Field
def B(x,y):
H11 = 4/((L+2*y)*np.sqrt((4*(d-2*x)**2/(L+2*y)**2+4)))

H12 = 4/((L+2*y)*np.sqrt((4*(d+2*x)**2/(L+2*y)**2+4)))

H13 = 4/((L-2*y)*np.sqrt((4*(d+2*x)**2/(L-2*y)**2+4)))

H14 = 4/((L-2*y)*np.sqrt((4*(d-2*x)**2/(L-2*y)**2+4)))

H21 = 1/((d-2*x)**2*(np.sqrt(4*(d-2*x)**2/(L-2*y)**2+4)))

H22 = 1/((d-2*x)**2*(np.sqrt(4*(d-2*x)**2/(L+2*y)**2+4)))

H23 = 1/((d+2*x)**2*(np.sqrt(4*(d+2*x)**2/(L-2*y)**2+4)))

H24 = 1/((d+2*x)**2*(np.sqrt(4*(d+2*x)**2/(L+2*y)**2+4)))

H1 = (H11 - H12 - H13 + H14)
H2 = (-4*(d-2*x))*(H21 - H22) + (-4*(d+2*x))*(H23 - H24)

B1 = (H1*mu_0)*2666
B2 = (H2*mu_0)*1700
B = np.sqrt(B1**2+B2**2)
return np.dstack((B1,B2,B))

B1 = B(x,y)[:,:,0]
B2 = B(x,y)[:,:,1]
Bt = B(x,y)[:,:,2]

#Contour Plot
levels = np.logspace(-2.5,1.23,10,base=10)
cp = plt.contourf(x, y, Bt, levels=levels, cmap=plt.cm.plasma,extend='max')
cb = plt.colorbar(cp)
plt.clim(0,2)

#X and Y increments
dxy = [0.023/151,0.02/151]

#Gradients of B, Bx & By

#Vector Field Plot
x, y = np.mgrid[-1.15:1.15:20j, -1:1:20j]*1E-2
Bq = B(x,y)[:,:,2]
dxy = [0.023/20, 0.02/20]

quiv = plt.quiver(x, y, dBq[0], dBq[1])
plt.show()

plt.figure()
x = np.linspace(-0.0115,0.0115,151)
plt.plot(x,B2[75,:])
plt.title('Magnetic Field vs. Position in x-direction')
plt.xlabel('Position (m)')
plt.ylabel('Magnetic Field (T)')
plt.show()

plt.figure()
y = np.linspace(-0.01,0.01,151)
plt.plot(y,-B1[:,75])
plt.title('Magnetic Field vs. Position in y-direction')
plt.xlabel('Position (m)')
plt.ylabel('Magnetic Field (T)')
plt.show()

#-----------------------------------------

yslice = 0.00375

v_x = np.zeros(N)
v_y = np.zeros(N)

#Arrays to store counts
#count1 = [0]
#count2 = [0]
#count3 = [0]
#count4 = [0]
#count5 = [0]
count6 = [0]

for i in range(0,N):
while xp[i]**2 + yp[i]**2 > radius**2:

for i in range(0,N):
if -yslice < yp[i] < yslice:
count6[0] = count6[0] + 1

fig1 = plt.figure()
plt.plot(xp,yp,".")

Fmag_x = np.zeros(N)
Fmag_y = np.zeros(N)
Fdrag_x = np.zeros(N)
Fdrag_y = np.zeros(N)
deltaV_x = np.zeros(N)
deltaV_y = np.zeros(N)

timestep = 1

#Method to animate simulation
def update(j,plt,fig1):
global xp, yp, v_x, v_y, t, timestep, count
plt.cla()
if min(xp**2 + yp**2) < radius**2:
count6.append(0)
t6.append(t6[timestep-1] + deltaT)
print(t6[timestep])
#print('fmag: ' + str(Fmag_x[0]))
#print('fdrag: ' + str(Fdrag_x[0]))
print('v(x): ' + str(v_x[0]))
for i in range(0,N):
if xp[i]**2 + yp[i]**2 < radius**2:

#Interpolation Method
Ix = int((xp[i]-(-0.0115))/(0.0115-(-0.0115))*150)
#print(Ix)
Iy = int((yp[i]-(-0.01))/(0.01-(-0.01))*150)
#print(Iy)
dBx1 = dBx[Ix,Iy]
dBx2 = dBx[Ix,Iy+1]
dBx3 = dBx[Ix+1,Iy]
dBx4 = dBx[Ix+1,Iy+1]
dBy1 = dBy[Ix,Iy]
dBy2 = dBy[Ix,Iy+1]
dBy3 = dBy[Ix+1,Iy]
dBy4 = dBy[Ix+1,Iy+1]
T = (xp[i]-x[Ix])/(x[Ix+1]-x[Ix])
#print(T)
U = (yp[i]-y[Iy])/(y[Iy+1]-y[Iy])
#print(U)
dBxx = (1-T)*(1-U)*dBx1+(1-T)*U*dBx2+T*(1-U)*dBx3+T*U*dBx4
#print('dBx ' +str(dBxx))
dByy = (1-T)*(1-U)*dBy1+(1-T)*U*dBy2+T*(1-U)*dBy3+T*U*dBy4
#print('dBy ' +str(dByy))
#Magnetic Force in X and Y
Fmag_x[i] = moment*dBxx*(B(xp[i],yp[i])[:,:,0]/np.abs(B(xp[i],yp[i])[:,:,2]))
Fmag_y[i] = moment*dByy*(B(xp[i],yp[i])[:,:,1]/np.abs(B(xp[i],yp[i])[:,:,2]))

#Drag Force in X and Y

#Increase in velocity
deltaV_x[i] = (Fmag_x[i] - Fdrag_x[i])*deltaT/m
deltaV_y[i] = (Fmag_y[i] - Fdrag_y[i])*deltaT/m

#Velocity
v_x[i] = v_x[i] + deltaV_x[i]
v_y[i] = v_y[i] + deltaV_y[i]

#New X and Y coordinates of particles
xp[i] = xp[i] + (v_x[i]*deltaT)
yp[i] = yp[i] + (v_y[i]*deltaT)

#Counter
if -yslice < yp[i] < yslice:
count6[timestep] = count6[timestep] + 1
timestep = timestep + 1
plt.plot(xp,yp,".")

ani = animation.FuncAnimation(fig1, update, frames=1, fargs=(plt, fig1), interval=0.1)

#Distribution of particles
#theta = np.arctan2(yp,xp)
#plt.hist(theta,50)
#plt.title('Distribution of Particles Generated')
#plt.ylabel('Number of Particles')
#plt.show()

• Because of the amount of code you've provided (and how little time I have to go over it), my suggestion would be to use Pythons time module (import time) and 'wrap' each segment of your code with it to time how long it takes to compute said segment i.e start = time.time(), your_code_segment_here, print(time.time()-start) and run the code to see what takes the most time, then see if you can optimise the segments. Also, the overhead from calling $B(x,y)$ many times (like $\mathcal{O}(10^5)$ calls) will significantly slow down the computation. – mattos Apr 27 '19 at 8:40
• Would you mind writing down the equations you are using? – nicoguaro Apr 27 '19 at 18:15
• A brief look on your code and thinking on the moving particles problem, perhaps the problem is related with close distance interactions (collisions). One way of trying to avoid problems could be to check if a collision is about to occur (distance between a set of particles smaller than a given amount) and then integrate in smaller time steps only on that set of particles. This will only be beneficial if collisions are rare events. – Ertxiem - reinstate Monica Apr 28 '19 at 0:08
• @nicoguaro, the equations for the drag force is simply Stoke's Law for viscous fluids, the magnetic field equation is quite long so I'll edit the post to include an image of it. But then the field gradient is found using the np.gradient function and interpolating between points, this then allows the magnetic force to be calculated using moment x dB (dBxx in code). Then by subtracting these forces from each other the overall force is found, from which the increase in velocit can be calculated by dividing by the mass and multiplying by the time step. Then the position is found using xf=xi+vt. – sgehugh2 Apr 28 '19 at 19:31
• Also @Ertxiem this simulation assumes no particle collisions. – sgehugh2 Apr 28 '19 at 19:31

The big derivatives may be a problem when integrating numerically. It seems to me that you're using Euler's method with a fixed time step. You may consider other methods such as the predictor–corrector methods that may help you to get smaller integration errors.

Furthermore, in your model, you get big derivatives whenever the particles are close, since the forces due to magnetic interactions decrease with the distance between particles. Hence, although you do not have physical collisions, whenever particles are close enough you have a strong interaction between them. And larger forces mean that the change in velocity is higher within each time step.

So, my advice is that you compute the total acceleration that each particle would feel (you already have to do that in your code) and to treat separately the subset of particles that have an acceleration above a given threshold. For that subset of particles you can subdivide the time step and make more intermediate computations regarding the movement of those particles.

You may also consider to do the same thing for the subset of particles with velocities above another threshold.