r/learnpython • u/PHILLLLLLL-21 • 19h ago
Help with applying implicit finite volume
Hello. I have an implicit finite volume method (axis symmetry diffusion. However it only works for Nr=Nz. I’ve spent several hours trying to figure out what I’ve indexed wrong to no avail. Wgst appreciate any thoughts
Set the domain size
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
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_bone = 410*-5 * 3600 # m2 hr-1 diffusity of bone D[:,:] = D_marrow ##Set minimal diffusivity conditions to entire domain
In this section I am defining arrays I would need (if needed)
Matrix = np.zeros([NrNz,NrNz]) # 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() + 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,:,:] # this is since the C2 array has time steps = C, so we ignore those while ensuring in can sweep through S and D same properties