crank3D

crank3D

"""Solve the 3D diffusion equation using CN and finite differences."""
from time import sleep

import numpy as np
import matplotlib.pyplot as plt
import networkx as nx
from pylab import *

# The total number of nodes
nodx = 3
nody = 3
nodz = 3


nnodes = nodx*nody*nodz

#this is for the plotting
xmin = 0.0
xmax = 1.0
ymin = 0.0
ymax = 1.0
# Ny = 4
# Nx = 4
# tmin = 0.0
# tmax = 1000.0
# Nt = 3000

# The total number of times
ntimes = 10000
# The time step
dt = 0.5

# The spatial mesh size setting it as a 1 meter square box
h = 1.0


tmin = 0
tmax = dt*ntimes

x, dx = np.linspace(xmin, xmax, nodx, retstep=True)
y, dy = np.linspace(ymin, ymax, nody, retstep=True)
t, dt = np.linspace(tmin, tmax, ntimes, retstep=True)

G = nx.grid_graph(dim=[nodx,nody,nodz])
L = np.matrix(nx.laplacian(G))

#making an expression for the rod in the center of the device, just need
#to vary the diffusion constant
# the diffusion constant for wet sand is about .6 and for aluminum is about 100
D = 1*np.matrix(np.eye(nnodes,nnodes))


#this part is a bit tricky, trying to insert a heat rod
#in the center of the device and figuring out which node
#of the graph corresponds to which spatial location
#takes a bit of leg work. the Spatial symmetry of
#the rectangular box makes it not too bad though
# the rod will be a little off center like this
# but the coding is already dense as it is and
# that shouldnt affect the actual experimental
#value of it.

w = nodx*nody/2
#rodthick is the number of nodes thick that we want
#the aluminum rod to be
rodthick = 2
d = 0

for i in range(nodz):
	
	# this loop determines how thick the rod will be horizontally
	for j in range(rodthick):
		d = j*nodx
		#this one how thick vertically
		for r in range(rodthick):
			
		#the targeted index k for the diffusion matrix
			k = w + r + d
			D[k,k] = 1
			
			k = w
		#just stepping up one graph row now so as to spatially go up 1 node 
		#in the physical vertical
		
		
		
		
	w = w + nodx*nodz



#making an expression for the heat source to go into the rhs section
C = np.matrix(np.zeros((nnodes,nnodes)))

C = np.matrix(np.zeros((nnodes,nnodes)))
C[nnodes/2,nnodes/2] = 0

	

# The rhs of the diffusion equation
rhs = -D*L/h**2 + C

# Setting initial temperature
T = 60*np.matrix(np.ones((nnodes,ntimes)))
for i in range(nnodes/2):
    T[i,0] = 0;

# Setup the time propagator. In this case the rhs is time-independent so we
# can do this once.
ident = np.matrix(np.eye(nnodes,nnodes))
pmat = ident+(dt/2.0)*rhs
mmat = ident-(dt/2.0)*rhs
propagator = np.linalg.inv(mmat)*pmat

# Propagate  E is for energy conservation
E = np.zeros(ntimes)
for i in range(ntimes-1):
	E[i] = sum(T[:,i])
	T[:,i+1] = propagator*T[:,i]

# To plot 1 time
print E[1]


#need to convert the big string T into a matrix for plotting and visualization


# a 4D array 
t = np.zeros((nodx, nody, nodz, ntimes))


w = 0
for p in range(ntimes):
	for i in range(nodx):
		for j in range(nody):
			for r in range(nodz):
				t[i,j,r,p] = T[w, p]
				w = w + 1
	w = 0		
	

# print t[:,:,2, 60]


V, dV = np.linspace(0, 70, 21, retstep=True)

print V
CS = plt.contourf(x,y,t[:,:,1,900], V)

plt.ylabel('distance (m)')
plt.xlabel('distance (m)')

cbar = colorbar(CS)
cbar.ax.set_ylabel('Temperature (K)')

plt.show()