1

I have seen other people asking something similar but I can not figure out the problem anyway. I am trying to translate a Matlab code to Python and I have a problem after the following line in the FOR loop:

dx = abs(np.diff(g_coord(num)))

Below you have the code up to that loop. Any help will be appreciated. I really tried to fix it by myself but unsuccessfully. Sorry if it is a stupid mistake. The MATLAB lines are kept as Python comments in case it helps.

import numpy as np
from scipy.sparse import lil_matrix 

# physical parameters
seconds_per_yr = 60*60*24*365; # number of seconds in one year
lx = 10000 ; #length of spatial domain (m)
Cp = 1e3 ; # rock heat capacity (J/kg/K)
rho = 2700 ; # rock density (kg/mˆ3)
K = 3.3 ; # bulk thermal conductivity (W/m/K)
kappa = K/(Cp*rho); # thermal diffusivity (mˆ2/s)
Tb = 0 ; # temperatures at boundaries (o C)
A = 2.6e-6 ; # heat production (W/mˆ3)
H = A/(rho*Cp); # heat source term (o K/s) % numerical parameters
dt = 1000*seconds_per_yr ; # time step (s)
ntime = 5000 ; # number of time steps
nels = 40 ; # total number of elements
nod = 2 ; # number of nodes per element
nn = nels+1 # total number of nodes
dx = lx/nels ; # element size

g_coord =  np.arange(0, lx+1, dx)#[0:dx:lx] 
bcdof = np.array([1, nn]); #[ 1 nn ] ; boundary nodes
bcval = np.array([Tb, Tb]); #[ Tb Tb ] ; # boudary values

g_num = np.zeros((nod, nels), float); #zeros(nod,nels) ;
g_num[0,:]=np.arange(1, nn); #g_num(1,:) = [1:nn-1] ;
g_num[1,:]=np.arange(2, nn+1); #g_num(2,:) = [2:nn] ;


# initialise matrices and vectors
ff = np.zeros((nn,1), float); # system load vector
b = np.zeros((nn,1), float); # system rhs vector
lhs=lil_matrix((nn, nn)) #lhs = sparse(nn,nn); system lhs matrix
rhs=lil_matrix((nn, nn)) #rhs = sparse(nn,nn); system rhs matrix
displ = np.zeros((nn,1), float); # initial temperature (o C)


#-----------------------------------------------------
# matrix assembly
#-----------------------------------------------------

# Matlab version of the loop
#-----------------------------------------------------
#for iel=1:nels # loop over all elements
#    num = g_num(:,iel) ; # retrieve equation number
#    dx = abs(diff(g_coord(num))) ; # length of element
#    MM = dx*[1/3 1/6 ; 1/6 1/3 ] ;# mass matrix
#    KM = [kappa/dx -kappa/dx ; -kappa/dx kappa/dx ]; #diffn matrix
#    F = dx*H*[1/2 ; 1/2] ; # load vector
#    lhs(num,num) = lhs(num,num) + MM/dt + KM ; # assemble lhs
#    rhs(num,num) = rhs(num,num) + MM/dt ; # assemble rhs
#    ff(num) = ff(num) + F ; # assemble load
#end # end of element loop

#Python version of the loop
#-----------------------------------------------------
for iel in range(0, nels): # loop over all elements
    num = g_num[:,iel]  # retrieve equation number
    #print(num)
    dx = abs(np.diff(g_coord[num]))  # length of element
    MM = dx*(np.array([[1/3, 1/6],[1/6, 1/3]]))  # mass matrix
    KM = np.array([[kappa/dx, -kappa/dx],[-kappa/dx, kappa/dx]]) 
    F = dx*H*(np.array([1/2, 1/2])).reshape(-1,1) # load vector
    lhs[num,num] = lhs[num,num] + MM/dt + KM # assemble lhs
    rhs[num,num] = rhs[num,num] + MM/dt # assemble rhs
    ff[num] = ff[num] + F # assemble load
6
  • If I modify that line for this "dx = abs(np.diff(g_coord[num]))" using square brackets this other error message appears: "arrays used as indices must be of integer (or boolean) type" Commented Jul 25, 2017 at 19:47
  • Please clarify. The callable error occurred when you tried to index an array with () (matlab style), right? You correct that, and you get this integer error? New versions of numpy have gotten picky about the index type. Why did you intialize g_num as float? Commented Jul 25, 2017 at 21:34
  • Regarding the lil sparse matrix - this is a start. But collecting element values in coo style inputs is probably faster. At least that's how I used to do with in MATLAB. Duplicates in coo inputs are summed when converted to csr for calculation purposes. Commented Jul 25, 2017 at 21:40
  • Matlab use () for both function arguments and indexing. Python uses [] for indexing. Commented Jul 25, 2017 at 22:00
  • yes, you are right, maybe I made a mistake initializing g_num as float as well. I used lil sparse matrix bc it looked similar to the sparse() common in Matlab, but I did not exactly knew the differences with coo style. Thanks for your help and time. Commented Jul 26, 2017 at 20:07

1 Answer 1

0

The error seems to be because num is a float. Simply do:

dx = abs(np.diff(g_coord[np.int32(num)]))

However it raises another error a few lines later because num is a 2-element array. You know what the code should do which I do not. If you have more issues, you can comment below or edit your question with the first problem solved.

Also, I noticed that you left all the ; at the end of the lines as in Matlab. You do not need this in Python. Also, I think there is no need for you to specify float when you create the matrices of zeros, they naturally are float.

Sign up to request clarification or add additional context in comments.

1 Comment

thank you very much, I really appreciate you took the time to look at it, best regards.

Your Answer

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

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.