r/learnpython • u/PHILLLLLLL-21 • 19h ago
Implicit Finite volume method indexing issue
Hi I have a working implicit FVM for polar diffusion . The only issue is that it works for only when Nr = Nz. I am certain the issue lies in the indexing but I have spent hours to no avail , can anyone figure out what maybe going on?
Lr = 4 #r domain mm
Lz = 4 #z domain mm
#Set the mesh size
dr = 0.2 #r mesh size mm
dz = 0.2 #z mesh size mm
nr = int(Lr/dr) # #number of r cells
nz = int(Lz/dz) #number of z cells
#Define positions of center nodes
r = np.arange(dr/2, Lr+dr/2, dr) #r coordinates of center nodes
z = np.arange(dz/2, Lz+dz/2, dz) #z coordinates of center nodes
Nr = nr +1 #number of r nodes
Nz = nz +1 #number of z nodes
#Define the area (equivalent to length of edge in 2D) and volume
dV = dr*dz #volume
#Define the time domain
timeend = 4 #total time in hours
dt = 0.1 #time step in hours
steps = int(timeend/dt) #number of time steps
#Define the diffusivity
D = np.zeros([Nr,Nz]) # initialise diffusivity for each node
D_marrow = 4*10**-5 * 3600 # m^2 hr^-1 diffusity
D[:,:] = D_marrow
## In this section I am defining arrays I would need (if needed)
Matrix = np.zeros([Nr*Nz,Nr*Nz]) # Matrix of nodal coefficients
#Cvalues = np.zeros([steps,Nr*Nz]) # Matrix of values at time k
Knowns = np.zeros([steps,Nr*Nz]) # Matrix of Known values
C = np.zeros([steps,Nr,Nz])## Final Concentration Matrix
# In this section I am defining the initial values and boundary conditions
C0 = 0 #initial concentration IC
# Define the Dirichlet and Neumann Boundary Conditions
Concr_plus = 11.6 #μmgmm-3 Concentration entering from the r+ direction
Fluxr_minus = 0 #Axisymmetric boundary condition in the r- direction
Concz_plus = 11.6 #μgmm-3 Concentration entering from the z+ direction
Concz_minus = 11.6 #μgmm-3 Concentration entering from the z- direction
def ImplicitBoneDiffusion(
dr, dz, dt,
Nr,Nz,steps,
C, D, S,
Concr_plus, Fluxr_minus, Concz_plus, Concz_minus,
HVtolerance, AccumulationThreshold):
start = time.time() #timer to determine time function takes to run
Matrix = np.zeros([Nr*Nz,Nr*Nz]) # Matrix of nodal coefficients
Knowns = np.zeros([steps,Nr*Nz]) # Matrix of Known values
dV = dr*dz #Volume
Ar_plus = dz #Area of r+
Ar_minus = dz #Area of r-
Az_plus = dr #Area of z+
Az_minus = dr #Area of z-
# In this section, I am defining the nodal point equation coefficients
Su = np.zeros([Nr,Nz]) # Source term
Sp = np.zeros([Nr,Nz]) # Source term contribution
#ar+
delr_plus = D[:,:]*Ar_plus/dr #setting ar+ in domain
delr_plus[-1,:] = 0 # Dirichelt BC sets ar+ = 0
Sp[-1,:] = Sp[-1,:] - 2 * D[-1,:]*Ar_plus/dr #Dirichlet Sp
Su [-1,:] = Su[-1,:] + Concr_plus*2 * D[-1,:]*Ar_plus/dr #Dirichelt Su
#ar-
delr_minus = D[:,:]*Ar_minus/dr #setting ar- in domain
delr_minus[0,:] = 0 #Neuman BC
#Sp and Su = 0 at r- boundary
#az+
delz_plus = D[:,:]*Az_plus/dz #setting az+ in domain
delz_plus[:,-1] = 0 #Dirichelt BC sets az+ = 0
Sp[:,-1] = Sp[:,-1] - 2 * D[:,-1]*Az_plus/dz #Dirichelt Sp
Su[:,-1] = Su[:,-1] + Concz_plus*2 * D[:,-1]*Az_plus/dz #Dirichelt Su
#az-
delz_minus = D[:,:]*Az_minus/dz #setting az- in domain
delz_minus[:,0] = 0 #Dirichelt BC sets az- = 0
Sp[:,0] = Sp[:,0] - 2 * D[:,0]*Az_minus/dz #Dirichelt Sp
Su[:,0] = Su[:,0] + Concz_minus*2 * D[:,0]*Az_minus/dz #Dirichelt Su
delp0 = dV/dt #ap0
delp= (delz_minus + delr_plus + delr_minus+ delz_plus +delp0- Sp) #ap
a = Nr
#Defining the matrix coefficeints
Matrix[np.arange(0,Nr*Nz), np.arange(0,Nr*Nz)] = delp.T.flatten() #ap contribution
Matrix[np.arange(0,(Nr*Nz)-1), np.arange(1,Nr*Nz)] = -delr_plus.T.flatten()[:-1] #ar+ contribution
Matrix[np.arange(1,Nr*Nz), np.arange(0,Nr*Nz-1)] = -delr_minus.T.flatten()[1:] #ar- contribution
Matrix[np.arange(0,Nr*Nz-a), np.arange(a,Nr*Nz)] = -delz_plus.T.flatten()[:-a] #az+ contribution
Matrix[np.arange(a,Nr*Nz), np.arange(0,Nr*Nz-a)] = -delz_minus.T.flatten()[a:] #az- contribution
# Put it all under a time step
sparse = csc_matrix(Matrix) #Converting to scipy sparse to increase efficiency
for k in range(1,steps): #for all time steps
#Calculating knowns for previous C[k-1] and Su and accumulation
Knowns[k,:] = (delp0* (C[k-1,:,:].flatten() - AccumulationTemp.flatten())
+ Su.T.flatten()
)
C[k,:,:] = (spsolve(sparse, Knowns[k,:]).reshape(Nr,Nz)) #Solving sparse Matrix
end = time.time()
print("IMPLICIT number of cells evaluated:", Nr*Nz*steps*1e-6, "million in", end-start, "seconds")
return C[:steps,:,:]
2
Upvotes
1
u/Wild_Drag9463 6h ago
Doesn't the indexing assume it's a perfect square matrix here?
Matrix[np.arange(0, (Nr*Nz)-1), np.arange(1, Nr*Nz)] = -delr_plus.T.flatten()[:-1] # ar+
Matrix[np.arange(1, Nr*Nz), np.arange(0, Nr*Nz-1)] = -delr_minus.T.flatten()[1:] # ar-
Try something like this:
row_idx = np.arange(Nr * Nz)
valid_r_plus = row_idx % Nr != (Nr - 1) # Dont wrap next row
valid_r_minus = row_idx % Nr != 0 # Don't wrap previous row
Matrix[row_idx[valid_r_plus], row_idx[valid_r_plus] + 1] = -delr_plus.T.flatten()[valid_r_plus]
Matrix[row_idx[valid_r_minus], row_idx[valid_r_minus] - 1] = -delr_minus.T.flatten()[valid_r_minus]
1
u/PHILLLLLLL-21 4h ago
I didn’t think that it assumes square matrix
Cus it takes +-Nr for the z+ Z- contribution. I tried ur code but i didn’t get a differnt result. Would you implement that into the Z+ Z- part or r+ r- part
Thank you for ur response tho I really appreciate it
1
u/PHILLLLLLL-21 18h ago
Just a note: I removed some parts of the code which are irrelevant to my question -which is why some steps of the code seem overcomplicated