#!/usr/bin/env python
# encoding: utf-8
"""
negative_mesh_forum_ex

Created by Yun Tao on 2012-04-01.
Copyright (c) 2012 __MyCompanyName__. All rights reserved.
"""

from fipy import *
from fipy import numerix
from matplotlib import pylab
import numpy as np
# mesh config
l = 50. 
nx = 800.
dx = l / nx
mesh = Grid1D(nx = nx, dx = dx) + [[-l/2]]


convection = VanLeerConvectionTerm 

# initial condition config
x = mesh.getCellCenters()[0]
distr = 1/numerix.sqrt(2*numerix.pi) * numerix.exp(-((x - 6) ** 2) / 2)

# nomralization of initial distribution
def_integ = np.trapz(distr.value, x)
print 'area under curve:', def_integ

# variable config
phi = CellVariable(name='$\phi$', mesh=mesh, value=distr)
print "initial cell volume: ", phi.getCellVolumeAverage() * l 
psi = CellVariable(name='$\psi$', mesh=mesh, value=0.1)

# Boundary Conditions
BCs = (FixedFlux(faces=mesh.getFacesRight(), value=0.),
	   FixedFlux(faces=mesh.getFacesLeft(), value=0.))

# viewer config
viewer = Viewer(vars=(phi, psi),
			    limits={'ymin': 0., 'ymax': 0.5},
			    title="$\phi$")


# equation config
D = 10.0
pos = FaceVariable(mesh=mesh, value=mesh.getFaceCenters(), rank=1)
c = (5.,) 
#c_switch = tuple(np.array((c)) * np.tanh(x))
#print cCoeff
mu = 0.3

#*np.tanh(40*pos))
eq1 = (TransientTerm(var=phi) == DiffusionTerm(coeff = D, var=phi) 
      + convection(coeff=(c / psi *np.tanh(40*pos)), var=phi) ) 

eq2 = (TransientTerm(var=psi) == ImplicitSourceTerm(var=phi) 
      + ImplicitSourceTerm(coeff=mu*phi, var=psi) )

eq = eq1 & eq2

# temporal config
dt = 0.001#1.0 * dx / abs(2)   # A = 1.0 seems to be the optimal config.
steps = 2000
#peak = [max(phi[:])]

# solve config
for steps in range(steps):
	print 'steps:', steps	
	eq.solve(boundaryConditions=BCs, dt=dt)
	print 'loc of min: ', mesh.getCellCenters()[:, numerix.argmin(phi)], '   log10: ', np.log10(min(phi))#, '\n'
	print 'loc of max: ', mesh.getCellCenters()[:, numerix.argmax(phi)], '   log10: ', np.log10(max(phi))#, '\n'
	
	print 'area under curve: ', np.trapz(phi.value, x)
	
	print 'cell volume: ', phi.getCellVolumeAverage() * l 
	print 'cell volume 2: ', phi.getCellVolumeAverage() * mesh.getCellVolumes().sum(), '\n'
#	peak.append(max(phi[:]))
#	print 'psi:', psi
#	print 'phi:', phi
	viewer.plot()
	
#print peak