Search code examples
pythonnumpyfluid-dynamics

perform derivative inside a loop


Hi I wrote the following code,

Z=1 #particle position
dz=0.1
time=1
dt=0.1
v=0
Fz=0
for time in np.arange(1, 10, dt):
#####simulation loop#######
    theta=np.arccos(-Z/R) #contact angle
    theta_e=((math.pi*110)/180) #equilibrium contact angle
    Z_e=-R*np.cos(theta_e)#equilibrium position of particle
    C=3.14*gamma*(R-Z_e) #additive constant


    Fsz= (gamma*math.pi*(Z-Z_e)**2)+(tau*2*math.pi*math.sqrt(R**2-Z**2))+C
    Fz=Fsz+(0.5*deltaF*np.sin((2*math.pi/lamda)*(Z-Z_e)-phi))#surface force
    w_a=gamma*lamda_m**2*(1-np.cos(theta_e)) #work of adhesion
    epsilon_z=2*math.pi*R*np.sin(theta)*mu*(nu/(lamda_m**3))*np.exp(w_a/KbT)#transitional drag
    epsilon_s=khi*mu*((4*math.pi**2*R**2)/math.sqrt(Ad))*(1-(Z/R)**2)
    epsilon=epsilon_z+epsilon_s
    Ft=math.sqrt(2*KbT*epsilon)*series #thermal force
    v=(-np.diff(Fz,Z)+Ft)/epsilon ##new velocity
    Z=v*dt #new position

I was trying to calculate dFz/dzbut it is gave me following error,

  File "C:/Users/mohammad.hossain1/Desktop/particle.py", line 62, in <module>
    v=(-np.diff(Fz,Z)+Ft)/epsilon ##new velocity

  File "C:\Users\mohammad.hossain1\AppData\Local\Continuum\anaconda3\lib\site-packages\numpy\lib\function_base.py", line 1924, in diff
    slice1[axis] = slice(1, None)

IndexError: list assignment index out of range

As my initial condition is Fz=0 Z=0and it changes with time I suppose to get dFz/dz. I have imported all necessary module and all variables are defined properly at the beginning of the code.But I got error while I have introduced the derivative. So most likely my approach is not going with argument. Is it possible to show me the mistake that I have done during my coding.


Solution

  • It's hard too say for sure even with the stack trace. Obviously numpy is trying to access a list called slice1 with an axis index that the list does not have. Not sure why this is the case, but it must be something with this line:

    v=(-np.diff(Fz,Z)+Ft)/epsilon
    

    I suspect specifically the part of this line that is causing this issue is the np.diff(), since it is np code throwing this error. My best guess is that Fz and Z are equal in this case, or otherwise unacceptable values for the .diff method. Try adding the following if like so:

    if Fz != Z:
         v=(-np.diff(Fz,Z)+Ft)/epsilon
    

    If that doesn't stop the crash try printing the values of Fz and Z right before this line and seeing if they look weird/suspicious.