-1
$\begingroup$

enter image description hereThis 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
radius = 0.0095
d = 0.023
L = 0.02

#Constants
mu_0 = 1.25663706E-6

#Particle Properties
p_rad = (4.9E-6)/2
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
dB = np.gradient(Bt,dxy[0],dxy[1])
dBy = np.gradient(B1,dxy[1])[0]
dBx = np.gradient(B2,dxy[0])[1]

#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]
dBq = np.gradient(Bq,dxy[0],dxy[1])

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

#X Gradient Plot
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()

#Y Gradient Plot
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

xp = np.random.uniform(-radius,radius,N)
yp = np.random.uniform(-radius,radius,N)

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:
        xp[i] = np.random.uniform(-radius,radius,1)
        yp[i] = np.random.uniform(-radius,radius,1)


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)
                #Interpolated Gradients
                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
                Fdrag_x[i] = 6*np.pi*mu*p_rad*v_x[i]
                Fdrag_y[i] = 6*np.pi*mu*p_rad*v_y[i]

                #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]
                #v_x = moment*30/(6*np.pi*mu*p_rad)*xp[i]/np.abs(xp[i])
                #v_y = moment*40/(6*np.pi*mu*p_rad)*yp[i]/np.abs(yp[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,".")
    plt.xlim(-radius,radius)
    plt.ylim(-radius,radius)

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.xlabel('Theta (Rads)')
#plt.ylabel('Number of Particles')
#plt.show()
Improve this question
$\endgroup$
5
  • $\begingroup$ 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. $\endgroup$
    – Matthew Cassell
    Commented Apr 27, 2019 at 8:40
  • $\begingroup$ Would you mind writing down the equations you are using? $\endgroup$
    – nicoguaro ♦
    Commented Apr 27, 2019 at 18:15
  • $\begingroup$ 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. $\endgroup$
    – Ertxiem - reinstate Monica
    Commented Apr 28, 2019 at 0:08
  • $\begingroup$ @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. $\endgroup$
    – sgehugh2
    Commented Apr 28, 2019 at 19:31
  • $\begingroup$ Also @Ertxiem this simulation assumes no particle collisions. $\endgroup$
    – sgehugh2
    Commented Apr 28, 2019 at 19:31

1 Answer 1

Reset to default
0
$\begingroup$

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.

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.