# -*- 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

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

#%% - INITIATION VALUES
C_a_0 = CellVariable(name = "initial value A", mesh=mesh) #Initial Value of A
C_b_0 = CellVariable(name = "initial value B", mesh=mesh) #Initial Value of B
C_c_0 = CellVariable(name = "initial value C", mesh=mesh) #Initial Value of C 
C_d_0 = CellVariable(name = "initial value D", mesh=mesh) #Initial Value of D
V0 = CellVariable(name = "initial velocity", mesh = mesh) #Initial Velocity 

C_a_0.setValue(0)
C_b_0.setValue(0)
C_c_0.setValue(0)
C_d_0.setValue(0)
V0.setValue(409.5183981003883) #V_flow/A_crossSection #m/hr

#%% - VARIABLES
C_a        = CellVariable(name="Concentration of A", mesh=mesh, hasOld=True, value = C_a_0)
C_b        = CellVariable(name="Concentration of B", mesh=mesh, hasOld=True, value = C_b_0)
C_c        = CellVariable(name="Concentration of C", mesh=mesh, hasOld=True, value = C_c_0)
C_d        = CellVariable(name="Concentration of D", mesh=mesh, hasOld=True, value = C_d_0)
V          = CellVariable(name="Velocity", mesh=mesh, hasOld=True, value=V0)


#%% - REACTION PARAMETERS
Da = 1 #Assumed diffusion Coefficient for all chemicals
k1 = 2000 #WAG

X = 0.8 #Conversion
Y = 0.8 #Yield

r_a = -k1*C_a

#%% - SET UP THE VIEWER
if __name__ == '__main__':
    viewer = Matplotlib1DViewer(vars=(C_a,C_b,C_c,C_d))
    viewer.plot()

#%% - BOUNDARY CONDITIONS
# Feeds (Left side BCs)
C_a_f = 3.19 #kmol/hr
C_b_f = (C_a_f*11)*1.15 #kmol/hr
C_c_f = 0 #kmol/hr
C_d_f = 0 #kmol/hr

C_a_BC= C_a_0*(1-X) #Assumed amount of A converted, reactor right BC

#FiPy Constraints
C_a.constrain(C_a_f, mesh.facesLeft)
C_a.constrain(C_a_BC,mesh.facesRight)
C_b.constrain(C_b_f,mesh.facesLeft)
C_c.constrain(C_c_f,mesh.facesLeft)
C_d.constrain(C_d_f,mesh.facesLeft)
V.constrain(V0,mesh.facesLeft)

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

EqC_a = (TransientTerm(var=C_a)==-ExponentialConvectionTerm(coeff=V_vec,var=C_a) + DiffusionTerm(var=C_b,coeff=Da) + r_a)

EqC_b = (TransientTerm(var=C_b)==-ExponentialConvectionTerm(coeff=V_vec,var=C_b) + DiffusionTerm(var=C_b,coeff=Da) + 11*r_a)

EqC_c = (TransientTerm(var=C_c)==-ExponentialConvectionTerm(coeff=V_vec,var=C_c) + DiffusionTerm(var=C_c,coeff=Da) + -17*r_a)

EqC_d = (TransientTerm(var=C_d)==-ExponentialConvectionTerm(coeff=V_vec,var=C_d) + DiffusionTerm(var=C_d,coeff=Da) + -8*r_a)

Eq = EqC_a & EqC_b & EqC_c & EqC_d


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

while elapsed < 10.:
    elapsed += dt
    C_a.updateOld()
    C_b.updateOld()
    C_c.updateOld()
    C_d.updateOld()
    res = 1.
    while res > 1.e-3:
        res = Eq.sweep(dt=dt, solver=solver)
    if __name__ == '__main__':
        viewer.plot()
        raw_input("Press <return> for next step...")