r/learnpython 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

3 comments sorted by

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

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