# -*- coding: utf-8 -*-
"""
Created on Wed Jul 13 14:05:25 2016

@author: ddesantis
"""

# -*- coding: utf-8 -*-
"""
Created on Tue Jul 12 13:32:09 2016

@author: ddesantis
"""
#%% - IMPORTS
import numpy as np
from fipy import *

#%% - REACTOR DIMENSIONS
d      = 0.4188 #m, diameter
L      = 1.2563 #m, length
r      = d / 2  #m, radius
A_crossSection = np.pi * r**2 #m^2, cross section
volFlow = 50 #m^3/hr

#%% - MESH
nx = 100
dx = L / nx
mesh = Grid1D(nx=nx, dx=dx)
x = mesh.cellCenters

#%% - INITIATION VALUES
C_a_0 = 0
C_b_0 = 0
C_c_0 = 0
C_d_0 = 0
C_d_0 = 0
C_e_0 = 0
T0 = 300+273.5 #K
V0 = volFlow/A_crossSection/60    #m/min


#%% - VARIABLES
C_a        = CellVariable(name="Concentration of A", mesh=mesh, hasOld=True)
C_b        = CellVariable(name="Concentration of B", mesh=mesh, hasOld=True)
C_c        = CellVariable(name="Concentration of C", mesh=mesh, hasOld=True)
C_d        = CellVariable(name="Concentration of D", mesh=mesh, hasOld=True)
C_e        = CellVariable(name="Concentration of E", mesh=mesh, hasOld=True)
V          = CellVariable(name="Velocity", mesh=mesh, hasOld=True)
T          = CellVariable(name="Temperature", mesh=mesh, hasOld=True)

C_a.setValue(C_a_0)
C_b.setValue(C_b_0)
C_c.setValue(C_c_0)
C_d.setValue(C_d_0)
C_e.setValue(C_e_0)
V.setValue(V0)
T.setValue(T0)

#%% - REACTION PARAMETERS
"""Reactions"""
# Simultaneous Reactions are occuring with compound a being the primary reagent.
# The reaction occurs over a catalyst bed
# Reaction 1: A + 11B --> 17C + 8D
# Reaction 1: A --> C + 5B

"""Coefficients"""
Da      = 100   # Assumed diffusion Coefficient for all chemicals
keff    = 100   # WAG
k1      = 100   # WAG
k2      = 100    # WAG

""" Conversion and Yield """
X = 0.8 # Conversion
Y = 0.8 # Yield

""" Heats of Formation """

H_a     = -110.4*1000 # kJ/kmol
H_b     = -241.8*1000 # kJ/kmol
H_c     = 0*1000      # kJ/kmol
H_d     = -393.8*1000 # kJ/kmol
H_e     = 20.3*1000   # kJ/kmol

""" Heats of Reaction """
dHrxn_1 = (C_c*H_c + C_d*H_d) - (C_a*H_a + C_b*H_b)
dHrxn_2 = (C_e*H_e + C_b*H_b) - C_a*H_a

""" Catalyst Info """
rho_cat = 2650 # kg/m^3
Cp_cat = 0.92  # kJ/kg*C
MW_cat = 184.4 # kg/kmol
Cp_cat = Cp_cat*MW_cat #kJ/kmol*C   
#%% - BOUNDARY CONDITIONS
# Feeds (Left side BCs)
C_a_f = 3.19/volFlow #kmol/hr
C_b_f = (C_a_f*11)*1.15 #kmol/hr, Feed excess B (15%)
C_c_f = 0 #kmol/hr
C_d_f = 0 #kmol/hr
C_e_f = 0 #kmol/hr

# Constrains
C_a.constrain(C_a_f, where=mesh.facesLeft)
C_b.constrain(C_b_f, where=mesh.facesLeft)
C_c.constrain(C_c_f, where=mesh.facesLeft)
C_d.constrain(C_d_f, where=mesh.facesLeft)
C_e.constrain(C_e_f, where=mesh.facesLeft)
V.constrain(V0, where=mesh.facesLeft)
T.constrain(T0, where=mesh.facesLeft)

#%% - DEFINE EQUATIONS
V_vec=V[None]

Eq_a = (TransientTerm(var=C_a)== -ExponentialConvectionTerm(coeff=V_vec,var=C_a) + DiffusionTerm(var=C_a,coeff=Da) - ImplicitSourceTerm(var=C_a,coeff=1*k1)- ImplicitSourceTerm(var=C_a,coeff=k2))

Eq_b = (TransientTerm(var=C_b)== -ExponentialConvectionTerm(coeff=V_vec,var=C_b) + DiffusionTerm(var=C_b,coeff=Da) - ImplicitSourceTerm(var=C_a,coeff=11*k1) - ImplicitSourceTerm(var=C_a,coeff=5*k2*(1-Y)))

Eq_c = (TransientTerm(var=C_c)== -ExponentialConvectionTerm(coeff=V_vec,var=C_c) + DiffusionTerm(var=C_c,coeff=Da) + ImplicitSourceTerm(var=C_a,coeff=17*k1))

Eq_d = (TransientTerm(var=C_d)== -ExponentialConvectionTerm(coeff=V_vec,var=C_d) + DiffusionTerm(var=C_d,coeff=Da) + ImplicitSourceTerm(var=C_a,coeff=8*k1))

Eq_e = (TransientTerm(var=C_e)== -ExponentialConvectionTerm(coeff=V_vec,var=C_e) + DiffusionTerm(var=C_e,coeff=Da) + ImplicitSourceTerm(var=C_a,coeff=k2*(1-Y)))

EqT  = (TransientTerm(var=T,coeff=rho_cat*Cp_cat) + ExponentialConvectionTerm(var=T,coeff=V_vec*rho_cat*Cp_cat) == DiffusionTerm(var=T,coeff = keff) + ImplicitSourceTerm(var=C_a, coeff=k1*dHrxn_1) - ImplicitSourceTerm(var=C_a, coeff=k2*dHrxn_2*(1-Y)))

Eq = Eq_a & Eq_b & Eq_c & Eq_d & Eq_e & EqT


#%% - SET UP THE VIEWER
if __name__ == '__main__':
    viewer1 = Viewer(vars=(C_a,C_b,C_c,C_d,C_c))
    viewer2 = Viewer(vars = T)
    viewer1.plot()
    viewer1.plot()

#%% - SOLVE EQUATIONS
elapsed = 0.
dt = 1e-1
solver = LinearLUSolver(tolerance=1e-10)

while elapsed < 10.:
    elapsed += dt
    print elapsed
    C_a.updateOld()
    C_b.updateOld()
    C_c.updateOld()
    C_d.updateOld()
    C_e.updateOld()
    T.updateOld()
    res = 1.
    while res > 1.e-2:
        res = Eq.sweep(dt=dt, solver=solver)
    if __name__ == '__main__':
        viewer1.plot()
        viewer2.plot()
        
print "[[Temp]]"
print T