<jupyter_start><jupyter_text>Higher order eigenfunctions<jupyter_code># Testing initial condition function
L1 = 100
L2 = 100
ds = 1
s1 = s2 = np.arange(0, 100, ds)
gamma = 1. 
beta = 1.
m = 2
n = 3
Y, X = np.meshgrid(s1, s2)
Z = CmnNeumann(0, X, Y, L1, L2, gamma, beta, m, n).T
h = plt.contourf(X,Y,Z)
plt.ylabel('First MT head position \n $s_1$ (nm)')
plt.xlabel('Second MT head position \n $s_2$ (nm)')
plt.colorbar(h)
plt.show()
# Testing solver initialization
pde_solver_mn = FPPassiveParaCNSolver('tests/CN_tests/para_tests/CNnm_tests/CN23/FP_passive_params.yaml')
grid = CmnNeumannInit(pde_solver_mn, m,n)
print(grid.todense())
s1 = np.asarray(pde_solver_mn.s1)
#print(s1)
s2 = np.asarray(pde_solver_mn.s2)
cs = plt.contourf(s2, s1, grid.todense())
plt.ylabel('First MT head position \n $s_1$ (nm)')
plt.xlabel('Second MT head position \n $s_2$ (nm)')
plt.colorbar(cs)
plt.show()
pde_solver_mn.Save()
# Read in h5py file
h5_data = h5py.File('tests/CN_tests/para_tests/CNnm_tests/CN23/FP_pass_para_CN.h5', 'r')
L1 = h5_data.attrs['L1']
L2 = h5_data.attrs['L2']
gamma = h5_data.attrs['gamma']
beta = h5_data.attrs['beta']
dt = h5_data.attrs['dt']
# Get time points
time = np.asarray(h5_data['time'])
nt = np.asarray(time.size)
err_arr = np.zeros(nt-1)
rel_err_arr = np.zeros(nt-1)
# Get solution and space points
xl_dist = np.asarray(h5_data['XL_data/XL_distr'])
s1 = np.asarray(h5_data['MT_data/s1'])
s2 = np.asarray(h5_data['MT_data/s2'])

# Loop over time and get error for each step

Y, X = np.meshgrid(s1, s2)
for i in range(nt-1):
    Z = CmnNeumann(time[i], X, Y, L1, L2, gamma, beta,m,n)
    err = xl_dist[:,:,i] - Z
    #err_arr[i] = dt*np.sum(err)
    err_arr[i] = np.sqrt(dt*np.sum(np.square(err)))

time_point = 10
err_grid = xl_dist[:,:,time_point] - CmnNeumann(time[time_point], X, Y, L1, L2, gamma, beta,m,n)
#print(err_grid)
#print(np.sqrt(dt*np.sum(np.square(err_grid))))
print("Point in time: ", time_point)
print("Summed error at time: ", np.sqrt(dt*np.sum(np.square(err_grid))))
fig, axarr = plt.subplots(2,2, figsize=(14,10))
cs00 =axarr[0,0].contourf(X,Y,xl_dist[:,:,time_point])
cs01 = axarr[0,1].contourf(X,Y,CmnNeumann(time[time_point], X, Y, L1, L2, gamma, beta, m ,n))
cs10 = axarr[1,0].contourf(X,Y,err_grid)
fig.colorbar(cs00, ax=axarr[0,0])
fig.colorbar(cs01, ax=axarr[0,1])
fig.colorbar(cs10, ax=axarr[1,0])
axarr[1,1].scatter(time[:-1], err_arr)

plt.tight_layout()
# Plot number of crosslinks over time
xl_num_arr = np.zeros(nt)
n_tot = np.sum(xl_dist[:,:,0])
print("Total crosslink number: ", n_tot)
for i in range(nt):
    xl_num_arr[i] = (np.sum(xl_dist[:,:,i])/n_tot)-1

plt.plot(time, xl_num_arr)
plt.ticklabel_format(axis='y', style='sci', scilimits=(-2,2))
plt.ylabel("Relative crosslink \n number change")
plt.xlabel("Time")
plt.tight_layout()
# Get the flux at the edge of tubules (Minus and plus ends of both microtubules)
# This is a first order accurate calculation whereas the solver should be second order accurate

deriv_ds = np.asarray([[-1./ds,1./ds]])
print(deriv_ds.shape)
## S1 minus end flux
# Rearrange matrix slice to so you can perform matrix operations accurately
S1m_flux = np.swapaxes(xl_dist[0:2,:, :],0,2)
S1m_flux = np.swapaxes(S1m_flux,1,2)
# Take first order derivative
S1m_flux = np.squeeze(np.matmul(deriv_ds, S1m_flux)) # First index is time, second index is other tubule location

## S1 Plus end flux
S1p_flux = np.swapaxes(xl_dist[-2:,:, :],0,2)
S1p_flux = np.swapaxes(S1p_flux,1,2)
S1p_flux = np.squeeze(np.matmul(deriv_ds, S1p_flux))

## S2 minus end flux
S2m_flux = np.swapaxes(xl_dist[:,0:2, :],0,2)
S2m_flux = np.swapaxes(S2m_flux,1,2)
S2m_flux = np.squeeze(np.matmul(S2m_flux, deriv_ds.T)) 

## S2 Plus end flux
S2p_flux = np.swapaxes(xl_dist[:,-2:, :],0,2)
S2p_flux = np.swapaxes(S2p_flux,1,2)
S2p_flux = np.squeeze(np.matmul(S2p_flux, deriv_ds.T))

# Graphing
fig, axarr = plt.subplots(2,2, figsize=(14,10))
csm = axarr[0,0].contourf(time, s2, S1m_flux[:-1,:].T)
fig.colorbar(csm, ax=axarr[0,0])
csp = axarr[0,1].contourf(time, s2, S1p_flux[:-1,:].T)
fig.colorbar(csp, ax=axarr[0,1])
axarr[1,0].plot(time, np.sum(S1m_flux[:-1,:],axis=1))
axarr[1,1].plot(time, np.sum(S1p_flux[:-1,:],axis=1))
axarr[0,0].set_title("First MT, \n minus-end flux")
axarr[0,1].set_title("First MT, \n plus-end flux")
plt.tight_layout()

fig2, axarr2 = plt.subplots(2,2, figsize=(14,10))
cs2m = axarr2[0,0].contourf(time, s1, S2m_flux[:-1,:].T)
fig2.colorbar(cs2m, ax=axarr2[0,0])
cs2p = axarr2[0,1].contourf(time, s1, S2p_flux[:-1,:].T)
fig2.colorbar(cs2p, ax=axarr2[0,1])
axarr2[1,0].plot(time, np.sum(S2m_flux[:-1,:],axis=1))
axarr2[1,1].plot(time, np.sum(S2p_flux[:-1,:],axis=1))
axarr2[0,0].set_title("Second MT, \n minus-end flux")
axarr2[0,1].set_title("Second MT, \n plus-end flux")
plt.tight_layout()<jupyter_output>(1, 2)
StarcoderdataJupyterStructuredCleanDedup