This is for a project that I am working on. I am trying to model an intrusive dike of magma cooling in some host or country rock within the earths crust. I am fairly new to coding. I did my best to convert this code form another coding language to python. I have a basic idea of what is going on. I know that I am trying to index something out of range, but am not sure where and or how to fix it. Any help I could get would be greatly appreciated! Thanks in advance.
import numpy as np
import matplotlib.pyplot as plt
#Physical parameters
L = 100 #Length of modeled domain [m]
Tmagma = 1200 #Temp. magma [°C]
Trock = 300 #Temp. of country rock [°C]
kappa = 1e-6 #Thermal diffusivity of rock [m^2/s]
W = 5 #Width of dike [m]
day = 3600*24 #seconds per day
dt = 1*day #Timestep
print(kappa)
print(day)
print(dt)
#Numerical parameters
nx = 201 #Number of gridpoints in x-direction
nt = 500 #Number of timesteps to compute
dx = L/(nx-1) #Spacing of grid
x = np.linspace(-50,50,100) #Grid
print(dx)
print(x)
#Setup initial temp. profile
T = np.ones(np.shape(x))*Trock
T[x>=-W/2] = 1200
T[x>=W/2] = 300
print(T)
time = 0
for n in range(0,nt): #Timestep loop
#compute new temp.
Tnew = np.zeros((1,nx))
print(Tnew)
for i in range(2,nx-1):
Tnew[i] = T[i]+kappa*dt*((T[i+1]-2*T[i]+T[i-1])/(dt**2)) # Set BC's
Tnew[1] = T[1]
Tnew[nx] = T[nx]
#update temp. and time
T = Tnew
time = time+dt
#Plot solution
plt.figure()
plt.plot(x,Tnew)
plt.xlabel('x [m]')
plt.ylabel('Temerpature [°C]')
plt.legend()
surf = ax.plot_surface(X, Y, cmap=cm.coolwarm, linewidth=0, antialiased=False)
fig.colorbar(surf, shrink=0.5, aspect=5)
plt.show()
IndexError Traceback (most recent call last)
<ipython-input-51-e80d6234a5b4> in <module>()
37 print(Tnew)
38 for i in range(2,nx-1):
---> 39 Tnew[i] = T[i]+kappa*dt*((T[i+1]-2*T[i]+T[i-1])/(dt**2))
40
41 # Set BC's
IndexError: index 2 is out of bounds for axis 0 with size
Tnew = np.zeros((1,nx))
This will gives you one array of nx element inside array (Two dimensional array) and you need to access it like Tnew[0][i]
just change it to
Tnew = np.zeros((nx))
check numpy.zeros DOC
NOTE: After resolving this error you are going to face is Index error at "T[i+1]" which is because when your 'i' reaches at last element you are not going to get 'i+1' element in 'T'.