#!/usr/bin/env python
# encoding: utf-8
"""
untitled.py

Created by Yun Tao on 2012-04-29.
Copyright (c) 2012 __MyCompanyName__. All rights reserved.
"""
import CreateMovie as movie
import matplotlib.pyplot as plt
from fipy import *
from matplotlib import pylab
from fipy.tools.numerix import *
import sys
import os
#matplotlib.use('TKAgg')

# pylab config
fig = pylab.figure(figsize=(24, 12))

ax1 = pylab.subplot((241))
ax2 = pylab.subplot((242))
ax3 = pylab.subplot((243))
ax4 = pylab.subplot((244))
ax5 = pylab.subplot((245))
ax6 = pylab.subplot((246))
ax7 = pylab.subplot((247))

# mesh config
nx = 40
ny = nx
dx = 1./nx
dy = dx
L = dx * nx
mesh = Grid2D(dx=dx, dy=dy, nx=nx, ny=ny)

# variable config
phi1 = CellVariable(mesh = mesh, value = 0.)
phi2 = CellVariable(mesh = mesh, value = 0.)
phi3 = CellVariable(mesh = mesh, value = 0.)
phi4 = CellVariable(mesh = mesh, value = 0.)
phi5 = CellVariable(mesh = mesh, value = 0.)
phi6 = CellVariable(mesh = mesh, value = 0.)
phi7 = CellVariable(mesh = mesh, value = 0.)

# constants config
epsilon = 0.01
b = ((1,), (1,))

# equations config
eq1 = TransientTerm() == -CentralDifferenceConvectionTerm(coeff=b) + DiffusionTerm(coeff=epsilon)
eq2 = TransientTerm() == -ExponentialConvectionTerm(coeff=b) + DiffusionTerm(coeff=epsilon)
eq3 = TransientTerm() == -HybridConvectionTerm(coeff=b) + DiffusionTerm(coeff=epsilon)
eq4 = TransientTerm() == -PowerLawConvectionTerm(coeff=b) + DiffusionTerm(coeff=epsilon)
eq5 = TransientTerm() == -UpwindConvectionTerm(coeff=b) + DiffusionTerm(coeff=epsilon)
eq6 = TransientTerm() == -ExplicitUpwindConvectionTerm(coeff=b) + DiffusionTerm(coeff=epsilon)
eq7 = TransientTerm() == -VanLeerConvectionTerm(coeff=b) + DiffusionTerm(coeff=epsilon)

# boundary conditions
x, y = mesh.getFaceCenters()

facesTopRight = ((mesh.getFacesRight() & (y > L / 2))
                | (mesh.getFacesTop() & (x > L / 2)))
facesBottomLeft = ((mesh.getFacesLeft() & (y < L / 2))
                    | (mesh.getFacesBottom() & (x < L / 2)))

BCs = (FixedValue(faces=facesTopRight, value=0.),
       FixedValue(faces=facesBottomLeft, value=1.))


# viewer config
viewer1 = MatplotlibViewer(vars=phi1, 
                           title="Central Difference Convection",
                           axes=ax1,  # axes defines the position of plots
                           legend=None)

viewer2 = MatplotlibViewer(vars=phi2, 
                           title="Exponential Convection",
                           axes=ax2,  # axes defines the position of plots
                           legend=None)

viewer3 = MatplotlibViewer(vars=phi3, 
                           title="Hybrid Convection",
                           axes=ax3,  # axes defines the position of plots
                           legend=None)

viewer4 = MatplotlibViewer(vars=phi4, 
                           title="Power Law Convection",
                           axes=ax4,  # axes defines the position of plots
                           legend=None)

viewer5 = MatplotlibViewer(vars=phi5, 
                           title="Upwind Convection",
                           axes=ax5,  # axes defines the position of plots
                           legend=None)

viewer6 = MatplotlibViewer(vars=phi6, 
                           title="Explicit Upwind Convection",
                           axes=ax6,  # axes defines the position of plots
                           legend=None)

viewer7 = MatplotlibViewer(vars=phi7, 
                           title="Van Leer Convection",
                           axes=ax7,  # axes defines the position of plots
                           legend=None)


viewer = MultiViewer(viewers=(viewer1, viewer2, viewer3, viewer4, viewer5, viewer6, viewer7))

# solver config
timeStepDuration = 0.9 * dx / int(numerix.sqrt(numerix.array(b[0])**2 + numerix.array(b[1])**2))
# ^ follows cfl condition: dt <= cfl * dx / velocity

steps = 50

for step in range(steps):
	
	eq7.solve(var=phi7,
			 boundaryConditions = BCs,
			 dt=timeStepDuration)
	viewer7.plot()
	
	
	eq6.solve(var=phi6,
			 boundaryConditions = BCs,
			 dt=timeStepDuration)
	viewer6.plot()
	
	eq5.solve(var=phi5,
			 boundaryConditions = BCs,
			 dt=timeStepDuration)
	viewer5.plot()
	
	eq4.solve(var=phi4,
			 boundaryConditions = BCs,
			 dt=timeStepDuration)
	viewer4.plot()
	
	eq3.solve(var=phi3,
			 boundaryConditions = BCs,
			 dt=timeStepDuration)
	viewer3.plot()

	eq2.solve(var=phi2,
			 boundaryConditions = BCs,
			 dt=timeStepDuration)
	viewer2.plot()
	
	eq1.solve(var=phi1,
			 boundaryConditions = BCs,
			 dt=timeStepDuration)
	viewer1.plot()
	
	fname = '_tmp%05d.png'%step
	plt.savefig(fname)
	
os.system("ffmpeg -r "+str(10)+" -b 1800 -i _tmp%05d.png 7convections.mp4")
os.system("rm _tmp*.png")