#
# A freely-propagating premixed hydrogen/air or hydrogen/O2 flame 
#
#
from Cantera import *
from Cantera.OneD import *
from Cantera.OneD.FreeFlame import FreeFlame
import numpy as np
from matplotlib import pyplot as pl
import sys


################################################################
# !!! GRID DEPENDENCE !!!
# Convergence is performed via refinement parameters (do not modify !)
# Grid convergence is ensured by modifying the nb of points in the initial grid
nb_grid = 10
# First use 10 and then test with 50 and 200 to ensure that flame speed is well converged
# !!!!!!!!!!!!!!!!!!!!!!!

#
# parameter values
#
p          =   1.0e5          # pressure
tin        =   300.0          # unburned gas temperature
rxnmech    =  'NAJM.cti'            # reaction mechanism fil
mix        =  'NAJM'              # gas mixture model
flag       =  'AIR'
#flag       =  'O2'

#
# Set the number of equivalence ratios
#
nb_phi = 1
vecphi = np.linspace(1,1.5,nb_phi)
T = []

C_atoms = 1.0
H_atoms = 4.0
wN2 = 2.8000000E+01
wO2 = 3.2000000E+01
wC  = 1.2000000E+01
wH  = 1.0000000E+00
wF  = C_atoms*wC + H_atoms*wH
s   = 4.0000000E+00  

if flag == 'O2':
	a = 0.0000100E+00
elif flag == 'AIR':
	a = 3.7600000E+00

for k in range(nb_phi):
	phi = vecphi[k]
        print 'Phi == '+str(phi)
	yf = (phi/s) / (phi/s+1+a*wN2/wO2)
	yo =    1    / (phi/s+1+a*wN2/wO2)
	yn = 1-yf-yo
	w  =    1    /(yf/wF+yo/wO2+yn/wN2)
	xo = (w/wO2) * yo
	xf = (w/wF ) * yf
	xn = (w/wN2) * yn

        comp = 'CH4:'+str(xf)+', O2:'+str(xo)+', N2:'+str(xn)  # premixed gas composition

# The solution domain is chosen to be 50 cm, and a point very near the
# downstream boundary is added to help with the zero-gradient boundary
# condition at this boundary.

	tol_ss = [1.0e-5, 1.0e-13]
	tol_ts = [1.0e-9, 1.0e-12]
# [rtol atol] for time stepping
	
	loglevel  = 1                       # amount of diagnostic output (0
	                                    # to 5)
					    
	refine_grid = 1                     # 1 to enable refinement, 0 to
	                                    # disable 				   
	
################ create the gas object ########################
#
# This object will be used to evaluate all thermodynamic, kinetic,
# and transport properties
#
	gas = IdealGasMix(rxnmech, mix)
	
# set its state to that of the unburned gas at the burner
	gas.set(T = tin, P = p, X = comp)
	
	print gas
	
	initial_grid = np.linspace(0, 0.01, nb_grid)
	f = FreeFlame(gas = gas, grid = initial_grid, tfix=900)
	
# set the properties at the inlet
	f.inlet.set(mole_fractions = comp, temperature = tin)
	
	f.set(tol = tol_ss, tol_time = tol_ts)
	f.setMaxJacAge(10, 10)
	f.set(energy = 'off')
	f.init()
	f.showSolution()
	
	f.solve(loglevel, refine_grid)
	
	f.setRefineCriteria(ratio = 5.0, slope = 0.05, curve = 0.005, prune = 0.0)
	#f.setRefineCriteria(ratio = 10.0, slope = 1, curve = 1)
	f.set(energy = 'on')
	f.solve(loglevel,refine_grid)
	
	outputname = 'P'+str(p/1e5)+'bar_T'+str(tin)+'K_phi'+str(phi)+'_'+flag+'_'+mix
	#f.save(outputname+'.xml')
	f.showSolution()
	
# Create reaction list header
        Lreac   = list(gas.fwdRatesOfProgress())
        for index, object in enumerate(Lreac):
           Lreac  [index] = "rrate_" +str(index+1)
        
        Lreac_r = list(gas.revRatesOfProgress())
        for index, object in enumerate(Lreac_r):
           Lreac_r[index] = "rrate_r_" +str(index+1)

# write the velocity, temperature, and mole fractions to a CSV file
	z = f.flame.grid()
	T = f.T()
	u = f.u()
	V = f.V()
	fcsv = open(outputname+'.csv','w')
	writeCSV(fcsv, ['#z (m)', 'u (m/s)', 'V (1/s)', 'T (K)', 'rho (kg/m3)']
	         + list(gas.speciesNames())
                 + Lreac
                 + Lreac_r )
	for n in range(f.flame.nPoints()):
	    f.setGasState(n)
	    writeCSV(fcsv, [z[n], u[n], V[n], T[n], gas.density()]
	             +list(gas.massFractions())
                     +list( 1000*gas.fwdRatesOfProgress())
                     +list(-1000*gas.revRatesOfProgress()) )
	fcsv.close()
	
	print 'solution saved to '+outputname+'.csv'
	print 'flamespeed = ',u[0],'m/s'
	#print 'Enthalpy_mole: ', 1.e-6*gas.enthalpy_mole()
	
	f.showStats()