'''
Created on 17 apr 2013

@author: barkar
'''
"""
Import necessary modules
"""
import sys
sys.path.append('U:\\python')
from fipy import *
from numpy import *
import numpy as np
from ctypes import *
from tq import *
from s4_proc import *
import matplotlib
import matplotlib.cm as cm
import matplotlib.mlab as mlab
import matplotlib.pyplot as plt
from time import sleep


temp=773.15
press=101325.0

nel=2
nph=1

tdtdb = create_string_buffer("user FECRsub-assess.TDB", TC_STRLEN_MAX)
mobtdb = create_string_buffer("MOBFE2", TC_STRLEN_MAX)

phases_class = tc_phases_strings * nph
phases = phases_class()
phases[0].phase = 'BCC'

elements_class = tc_elements_strings * nel
elements = elements_class()
elements[0].element = 'CR'
elements[1].element = 'FE'
nval=s4_init(tdtdb,mobtdb,nel,elements,nph,phases,temp,press)

"""
Define system dimensions
"""
L = 1.0e-8
nx=250
dx=L/nx
mesh=Grid1D(nx=nx,dx=dx)
dist=np.linspace(0,L,nx+1)
xxx=np.linspace(0,0,nx+1)
yyy=np.linspace(0,0,nx+1)
x0=np.linspace(0,0,nx+1)

"""
Create a cellvariable and set initial values
"""
phi = CellVariable(name=r"$\phi$", mesh=mesh)
noice=GaussianNoiseVariable(mesh=mesh,mean=0.0,variance=0.0001)
noice2=GaussianNoiseVariable(mesh=mesh,mean=0.36,variance=0.00001)
PHI = phi.getArithmeticFaceValue()
x=mesh.cellCenters()


"""
Set initial noice
"""
phi.setValue(noice2)


"""
L from assessment
"""
L0 = 24212.06-15.507*temp
L1=1664.69+0.286*temp
a =2.5e-10
kappa=(L0*a**2)/2
R=8.3144
Vm=7.09*10e-6

     
"""
Create various vectors for storing thermodynamic data
"""
dmucrdcr=np.linspace(0,0,size(PHI))
dmucrdfe=np.linspace(0,0,size(PHI))
dmufedcr=np.linspace(0,0,size(PHI))
dmufedfe=np.linspace(0,0,size(PHI))
d2G=np.linspace(0,0,size(PHI))
M_Fe=np.linspace(0,0,size(PHI))
M_Cr=np.linspace(0,0,size(PHI))
M_Fe10=np.linspace(0,0,size(PHI))
M_Cr10=np.linspace(0,0,size(PHI))
  
L_CrFe=((PHI*(1-PHI)*(PHI*M_Fe+(1-PHI)*M_Cr)))


"""
Creating a viewer to be called when solvning the equations
"""
viewer = Viewer(vars=(phi,), datamin=0.0, datamax=1.0,colobar=True)

"""
Setting the time step, duration time, and solving the equations with updated thermodynamic
data for every iteration
"""

elapsed = 0.
duration = 100.0*3600
dexp = -5
dt=1.0
imgcounter=1

coeff1 = Variable(value=L_CrFe*d2G)
coeff21 = Variable(value=L_CrFe)
coeff22 = Variable(value=kappa)


eq = (TransientTerm() == DiffusionTerm(coeff=coeff1) - DiffusionTerm(coeff=(coeff21,coeff22)))

while elapsed < duration:

    elapsed += dt
    dexp += 0.5

    for i in range(size(PHI)):

        xxx[i]=PHI[i]
        valarr=s4_calc(nval,nel,xxx[i],temp,press)
        dmucrdcr[i]=valarr[2].value
        dmucrdfe[i]=valarr[3].value
        dmufedcr[i]=valarr[4].value
        dmufedfe[i]=valarr[5].value
        M_Cr[i]=pow(10,valarr[6].value)
        M_Fe[i]=pow(10,valarr[7].value)
        M_Cr10[i]=valarr[6].value
        M_Fe10[i]=valarr[7].value
        d2G[i]=dmucrdcr[i]-dmufedcr[i]    


    coeff1.value = L_CrFe*d2G
    coeff21.value = L_CrFe
    coeff22.value = kappa 


    
    for i in range(size(PHI)):
        yyy[i]=PHI[i]

    eq.solve(phi, dt=dt)
 
    for i in range(size(PHI)):
        xxx[i]=PHI[i]
        
    dd=0.0
    ddeps=0.000001

    for i in range(size(PHI)):
        dd=max(dd,abs(xxx[i]-yyy[i]))

    if (abs(dd-ddeps)<ddeps/10):
        dt=dt
    elif (dd < 0.001):
        dt=dt*1.5
    else:
        dt=dt/1.5

    viewer.plot()

 
    
    hours=round(elapsed/3600.)
    line="{0:7d} hours".format(int(hours))
    if ( hours < 1.0 ):
        minutes=round(elapsed/60.)
        line="{0:7d} minutes".format(int(minutes))

    plt.title(line, horizontalalignment='right')
    
    time=round(elapsed/3600)

    print "time:",time,"h    (dt :",dt,"s ",dd*100,"% change)"
            
    if elapsed>=duration:
        print "end of simulation, duration:",round(duration/3600./24.),"days (",round(duration/3600./24./365.),"years)"