currently I am using this code to output a 2-D implicit convection diffusion matrix. I cannot figure out how to get the last row of the column to populate. I have it grab all this data from a txt file which has also been included below the code. I then try to use the range values to have it out put a matrix but I cannot seem to get the matrix to output correctly. any help or an explanation as to why it wont work would be great appreciated. when I change the range of i and j to just Nx or Ny I get an error for indexing which I don't quite understand either, however in order to solve this have set the range to Nx-1 and Ny-1 which I believe to be incorrect. any advice is once again greatly appreciated.
import os
import numpy as np
#import data into variables for use
file = open("diffsolverin.txt","r")
line = file.readline() # rea the first line (remove comments from input)
tokens = line.split() # this splits the line into chunks of characters separated by spaces
Nx = int(tokens[0]) # assign the value of the first chunk of characters to Nx after interpreting them as an integer
Ny = int(tokens[1])
width = float(tokens[2])
height = float(tokens[3])
D = float(tokens[4])
t = float(tokens[5])
U_x = float(tokens[6])
U_y = float(tokens[7])
S_p = float(tokens[8])
line = file.readline() # read the second line
tokens = line.split()
left_bc_type = int(tokens[0]) #0 is Dirichlet, 1 is Neumann
right_bc_type = int(tokens[1])
top_bc_type = int(tokens[2])
bottom_bc_type = int(tokens[3])
line = file.readline() # read the third line
tokens = line.split()
left_bc_value = float(tokens[0])
right_bc_value = float(tokens[1])
top_bc_value = float(tokens[2])
bottom_bc_value = float(tokens[3])
N=Nx*Ny
Amat = np.zeros((N,N))
dx = Nx/width
dy = Ny/height
B = np.zeros(N)
phi = np.zeros(N)
#solve and establish boundary conditions
for i in range(Nx-1):
for j in range(Ny-1):
k = (j)*(Ny) i
print (k)
if i == 0: # this is a left boundary cell, deal with it given the BCs)
if left_bc_type == 0:
A_e = D*dx-U_x/2
A_w = 0
A_n = D*dy U_y/2
A_s = D*dy-U_y/2
A_p = A_e A_w A_s A_n
Amat[k,k] = A_p
Amat[k,k-Nx]=A_n
Amat[k,k Nx]=A_s
Amat[k,k 1]=A_e
B[k] = A_p*left_bc_value
elif left_bc_type == 1:
A_e = D*dx-U_x/2
A_w = 0
A_n = D*dy U_y/2
A_s = D*dy-U_y/2
A_p = A_e A_w A_s A_n
Amat[k,k] = A_p
Amat[k,k-Nx]=A_n
Amat[k,k Nx]=A_s
Amat[k,k 1]=A_e
B[k] = A_p
elif i == Nx-1: #this is a right boundary cell, ...
if right_bc_type == 0:
A_e = 0
A_w = D*dx U_x/2
A_n = D*dy U_y/2
A_s = D*dy-U_y/2
A_p = A_e A_w A_s A_n
Amat[k,k] = A_p
Amat[k,k-Nx]=A_n
Amat[k,k Nx]=A_s
Amat[k,k-1]=A_w
B[k] = A_p*right_bc_value
elif right_bc_type == 1:
A_e = 0
A_w = D*dx U_x/2
A_n = D*dy U_y/2
A_s = D*dy-U_y/2
A_p = A_e A_w A_s A_n
Amat[k,k] = A_p
Amat[k,k-Nx]=A_n
Amat[k,k Nx]=A_s
Amat[k,k-1]=A_w
B[k] = A_p
elif j == 0: # this is a top boundary cell
if top_bc_type == 0:
A_e = D*dx-U_x/2
A_w = D*dy U_x/2
A_n = 0
A_s = D*dy-U_y/2
A_p = A_e A_w A_s A_n
Amat[k,k] = A_p
Amat[k,k 1]=A_e
Amat[k,k Nx]=A_s
Amat[k,k-1]=A_w
B[k] = A_p*top_bc_value
elif top_bc_type == 1:
A_e = D*dx-U_x/2
A_w = D*dx U_x/2
A_n = 0
A_s = D*dy-U_y/2
A_p = A_e A_w A_s A_n
Amat[k,k] = A_p
Amat[k,k-1]=A_w
Amat[k,k Nx]=A_s
Amat[k,k 1]=A_e
B[k] = A_p
elif j == Ny-1: # this is a bottom boundary cell
if bottom_bc_type == 0:
A_e = D*dx-U_x/2
A_w = D*dx U_x/2
A_n = D*dy U_y/2
A_s = 0
A_p = A_e A_w A_s A_n
Amat[k,k] = A_p
Amat[k,k 1]=A_e
Amat[k,k-Nx]=A_n
Amat[k,k-1]=A_w
B[k] = A_p*bottom_bc_value
elif bottom_bc_type == 1:
A_e = D*dx-U_x/2
A_w = D*dx U_x/2
A_n = D*dy U_y/2
A_s = 0
A_p = A_e A_w A_s A_n
Amat[k,k] = A_p
Amat[k,k 1]=A_e
Amat[k,k-Nx]=A_n
Amat[k,k-1]=A_w
B[k] = A_p
else:
A_e = D*dx-U_x/2
A_w = D*dx U_x/2
A_n = D*dy U_y/2
A_s = D*dy-U_y/2
A_p = A_e A_w A_s A_n
Amat[k,k] = A_p
Amat[k,k 1]=A_e
Amat[k,k-Nx]=A_n
Amat[k,k-1]=A_w
Amat[k,k Nx]=A_s
B[k]=A_p
print (Amat)
#print(B)
#print (ans)
file = open("diffsolverout.txt", "w ")
filew = str(Amat)
filew1=str(B)
file.write(filew)
file.write(filew1)
file.close()
and this is the contents of the text file:
10 10 1.0 2.0 0.5 0.1 1 4 1
1 1 1 1
10 10 10 10
Once again I either get a blank last row of the matrix which I believe is incorrect or I get an index error for attempting to fill the matrix.
CodePudding user response：
You are reaching outside the boundaries of the matrix, as the error suggests. It is not from the value of i or j being wrong though. Your boundary cell definitions take k Nx or k - Nx, but when you inspect this it is wrong simply because of how indexing works for python.
The error occurs when k = 90, and you have Nx = 10, so k Nx = 90 10 = 100. But the indexing goes from 0-99 not 1-100. You should just be able to change each of these to be k (Nx - 1) or k - (Nx - 1) and fix this specific issue.
You were already able to identify this when saying that j == Ny - 1, it is the same logic.
CodePudding user response：
Ok so after looking into this more, the way you are indexing the Amat is incorrect for what you are trying to do. You have to understand how indexing in a matrix in python (and more generally in maths) works. I thought it deserved a new solution as I wanted to explain it properly.
There are i number of rows and j number of columns. Your approach is to label each cell in the matrix from left to right, top to bottom with an integer k. This is fine and helps out a lot in this approach, but you cannot outright use it to change values in the matrix the way you are. (As a side note, i is the rows and j is the columns, your approach, while setup to be this way, actually takes i as the columns and j as the rows when you define k)
Imagine that we have k = 5. Then in the matrix, we have A[5,5] as the point we want to deal with as the central point. To go left or right we get A[k, k-1] and A[k, k 1], but to go up or down we have to do A[k-1,k] and A[k 1,k].
For your else block (to show all examples) this should actually be:
Amat[k, k] = A_p
Amat[k, k 1] = A_e
Amat[k - 1, k] = A_n
Amat[k, k - 1] = A_w
Amat[k 1, k] = A_s
You'll notice that this just takes k as the point of reference but we still need to use the standard way of indexing to reach the values you want, not just add or subtract Nx to the cell number. If that makes sense.
Lastly, I just wanted to add a few suggestions to make your code more readable and a lot more debuggable.
The values for A_e, A_w, A_n, A_s can be defined just once at the start and then A_p can be defined using the required values for the specific case.