Crank3D

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

import numpy as np
import matplotlib.pyplot as plt
import networkx as nx

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

nnodes = nodx*nody*nodz
# The total number of times
ntimes = 500
# The time step
dt = 0.5
# The diffusion constant
D = 0.1
# The spatial mesh size
h = 1.0

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

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

# 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[2]
print T[:,100]

# To plot all times
#plt.plot(T)