#
# A freely-propagating premixed CH4/Oxygen 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 = 5
# 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    =  'FRASSO.cti'            # reaction mechanism fil
mix        =  'JL_R'              # gas mixture model
flag       =  'O2'

#
# Set the number of equivalence ratios
#
nb_phi = 1
vecphi = np.linspace(1.0,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)
	
	initial_grid = np.linspace(0, 0.01, nb_grid)
	f = FreeFlame(gas = gas, grid = initial_grid, tfix=1000)
	
# 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.restore(file = 'save.xml', id = '1')
#        f.showSolution()
#        sys.exit()
	f.solve(loglevel, refine_grid)
	
	f.setRefineCriteria(ratio = 400.0, slope = 0.4, curve = 0.04, prune = 0.0)
#	f.setRefineCriteria(ratio = 1.0, slope = 0.1, curve = 0.05, prune = 0.0)
	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()
	
	f.save(file='save.xml', id='1', desc='none')

# write the velocity, temperature, and mole fractions to a CSV file
	z = f.flame.grid()
	T = f.T()
	u = f.u()
	V = f.V()
# 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)
        
        Lsourcespec = list(gas.speciesNames())
        for index, object in enumerate(Lsourcespec):
           Lsourcespec[index] = "specsource_" +object
        
# Computes chem source terms in kg/m3/s
        nuki=gas.productStoichCoeffs()-gas.reactantStoichCoeffs()
        nspec=len( gas.speciesNames() )
        nreac=len( gas.netRatesOfProgress() )
        specsource = zeros((nspec,f.flame.nPoints()),'d')
        specsource_list = [ ]
        for n in range(f.flame.nPoints()):
          f.setGasState(n)
          for k in range(nspec):
            rrate_progress=1000*gas.netRatesOfProgress()
            #for i in range(nreac):
              # molW[k] is in g/mol
              #specsource[k,n]=specsource[k,n]+nuki[k,i]*molW[k]/1e3*rrate_progress[i]
              #print "nuki[k,i],rrate_progress[i]",nuki[k,i],rrate_progress[i]
              
          
          #specsource_list.append( [specsource[k,n] for k in range(nspec)] )
          #print "node ",n
          #print "specsource_list[n] ",specsource_list[n]
        specsource_list = numpy.array(specsource_list)

	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'
	
	f.showStats()