<jupyter_start><jupyter_text>This code compares net reaction rate computation of TChem++ with cantera<jupyter_code>def runTchem():
    os.chdir(dirr) 
    os.system(direx+"tchem-example-reactionrates-surface.x")
    os.chdir('cantera')
    return

def printRoP():
    dirOputs = dirInputs
    species_names =  gas.species_names + surf.species_names
    ## gas 
    omega = np.loadtxt(dirOputs+'omega.dat')
    #surf 
    omegaGasSurf = np.loadtxt(dirOputs+'omegaGasSurf.dat')
    omegaSurf    = np.loadtxt(dirOputs+'omegaSurf.dat')
    omegaS = np.vstack((omegaGasSurf,omegaSurf))
    print('--Surface phase --')
    for i,species in enumerate(species_names):
        print(species,surf.net_production_rates[i],omegaS[i,1])
        
    print('--Surface rate of progress--')
    for i,reac in enumerate(surf.reaction_equations()):
        print(reac,surf.forward_rates_of_progress[i],surf.reverse_rates_of_progress[i])
        
    print('--- reaction constatns ---')
    for i,reac in enumerate(surf.reaction_equations()):
        print(reac,surf.reverse_rate_constants[i],surf.forward_rate_constants[i])
    print('--Gas phase --') 
    for i,species in enumerate(gas.species_names):
        print(species,gas.net_production_rates[i],omega[i,1])
   
    return 


def compareTchemAndCantera():
    dirOputs = dirInputs
    net_prodGasPhase  = np.loadtxt(dirOputs+'omega.dat') [:,1] 
    net_prodGasSurf   = np.loadtxt(dirOputs+'omegaGasSurf.dat')[:,1] 
    net_prodSurfSurf  = np.loadtxt(dirOputs+'omegaSurf.dat')[:,1] 
    net_prodSurf = np.hstack((net_prodGasSurf,net_prodSurfSurf))
    
    print('-----------------------------')
    print('Gas phase')
    print('species name cantera  Tchem diff')
    print('-----------------------------')
    for i,species in enumerate(gas.species_names):
        print(species,gas.net_production_rates[i],net_prodGasPhase[i],
              gas.net_production_rates[i]-net_prodGasPhase[i],
             (gas.net_production_rates[i]-net_prodGasPhase[i])/(gas.net_production_rates[i]+1e-23))
    
        if abs((gas.net_production_rates[i]-net_prodGasPhase[i])/gas.net_production_rates[i]) > allow_error:
            print('CHECK this rate')
        print('    ')
        
    print('-----------------------------')
    print('Surface interface') 
    print('species, name, cantera,  Tchem,  diff, % diff')
    print('-----------------------------')
    
    species_names =  gas.species_names + surf.species_names
    for i,species in enumerate(species_names):
        print(species,surf.net_production_rates[i],
              net_prodSurf[i],surf.net_production_rates[i]-net_prodSurf[i],
             (surf.net_production_rates[i]-net_prodSurf[i])/(surf.net_production_rates[i]+1e-23))
    
        if abs((surf.net_production_rates[i]-net_prodSurf[i])/surf.net_production_rates[i]) > allow_error:
            print('CHECK this rate')
    
        print('    ')
    return 

 
def getOriginalInput():
    #Reading mass-fraction:
    #======================
    state = np.loadtxt(dirInputs+'inputGasOriginal.dat') # ρ,P,temp, mass fraction
                                  
    Ykk = np.hstack((np.atleast_1d(state[2]),np.atleast_1d(state[1]),state[3:]))
    #Reading site-fraction:
    #======================
    Zkk = np.loadtxt(dirInputs+'inputSurfGasOriginal.dat')
    return Ykk , Zkk

def saveInputsTChem():
    state =  np.hstack((np.atleast_1d(gas.density),np.atleast_1d(Ykk[1]),np.atleast_1d(Ykk[0]),Ykk[2:]))
    # T, P
    np.savetxt(dirInputs + 'inputGas.dat',state)
    np.savetxt(dirInputs + 'inputSurfGas.dat',Zkk)
    return 

def getRandomSample(Nvars,μ=0,σ=1,sum1=True):
    sample = np.random.rand(Nvars)*σ+μ
    if sum1:
        sample /= np.sum(sample)
    
    return sample

def makeTChemRandomSample():
    massFraction = getRandomSample(gas.n_species)
    T            = getRandomSample(1,μ=1200,σ=500,sum1=False)
    P            = 1013250.0 # Pascal
    Ykk = np.hstack((T,P,massFraction))
    Zkk = getRandomSample(surf.n_species)  
    return Ykk,Zkk<jupyter_output><empty_output><jupyter_text>User inputs* dirInputs : input and/or ouput * direx: Where TChem executable (tchem-example-reactionrates-surface.x) is located. source code is at /src/example/TChem_ReactionRatesSurface.cpp* dirr : run directory     * allow_error: allow error between cantera and Tchem* make_random_sample: true or false, if true a random state vector will be used. else a vector will be loaded inputGasOriginal.dat and inputSurfGasOriginal.dat.<jupyter_code>TchemRepo = '/Users/odiazib/CODE/TChem++/'
dirr      = TchemRepo+'build/example/SurfVerification/'
direx     = TchemRepo+'build/example/'
allow_error        = 1e-3
make_random_sample = True<jupyter_output><empty_output><jupyter_text>With an initial guess<jupyter_code># dirInputs = TchemRepo + 'build/example/SurfVerification/data/reaction-rates-surfaces/PT/'
# gas  = ct.Solution(dirInputs+'chem.cti', 'gas')
# surf = ct.Interface(dirInputs+'chem.cti', 'surface1', [gas])


dirInputs = TchemRepo + 'build/example/SurfVerification/data/reaction-rates-surfaces/X/'
gas  = ct.Solution(dirInputs+'chem_annotated-gas_surface.cti', 'gas')
surf = ct.Interface(dirInputs+'chem_annotated-gas_surface.cti', 'surface1', [gas])

if make_random_sample:
    Ykk , Zkk = makeTChemRandomSample()
else:    
    Ykk , Zkk = getOriginalInput()
#     Ykk , Zkk = getOriginalInput()

gas.Y  = Ykk[2:]
gas.TP = Ykk[0], Ykk[1]

surf.TP =Ykk[0], Ykk[1]
surf.coverages=Zkk
print('\nsurf.coverages after\n', surf.coverages)
saveInputsTChem()
runTchem() # run Tchem
print('*-------------*')
print('With an random guess')
compareTchemAndCantera()
# printRoP()<jupyter_output><empty_output><jupyter_text>Change gas temperature<jupyter_code>Ykk[0] = 2000

gas.Y  = Ykk[2:24]
gas.TP = Ykk[0], Ykk[1]
surf.TP = Ykk[0], Ykk[1]

saveInputsTChem()

runTchem()
print('*-------------*')
print('Changing gas temperature to '+str(Ykk[0])+' [K]')
compareTchemAndCantera()<jupyter_output>*-------------*
Changing gas temperature to 2000.0 [K]
-----------------------------
Gas phase
species name cantera  Tchem diff
-----------------------------
Ar 0.0 0.0 0.0 0.0
    
Ne 0.0 0.0 0.0 0.0
    
N2 0.0 0.0 0.0 0.0
    
CH4(2) 18747919.214019556 18748080.0 -160.7859804444015 -8.576204036774756e-06
    
O2(3) 135145334.00254986 135149700.0 -4365.997450143099 -3.230594294924546e-05
    
CO2(4) 0.0 0.0 0.0 0.0
    
H2O(5) 100052.57596896253 100056.4 -3.8240310374676483 -3.822021572591901e-05
    
H2(6) 2486971.3995598955 2487005.0 -33.60044010449201 -1.3510585650658499e-05
    
CO(7) 1410170.6742252484 1410198.0 -27.325774751603603 -1.9377636516669485e-05
    
C2H6(8) 2731845.2151357657 2731871.0 -25.78486423427239 -9.438625618835051e-06
    
CH2O(9) -4182.57296382031 -4182.92 0.3470361796898942 -8.29719368177897e-05
    
CH3(10) 116003763.6952291 116007900.0 -4136.304770901799 -3.565664284625201e-05<jupyter_text>Changing coverages<jupyter_code>if make_random_sample:
    Ykk , Zkk = makeTChemRandomSample()
else:    
    Ykk , Zkk = getOriginalInput()

Zkk[-1] *=2
Zkk[-2] *=0.5

gas.Y  = Ykk[2:24]
gas.TP = Ykk[0], Ykk[1]
#surf.Y=Zkk[2:19]
surf.TP = Ykk[0], Ykk[1] 
surf.coverages = Zkk 

Zkk /= np.sum(Zkk)
print(Zkk)   

print('\nsurf.coverages after\n', surf.coverages)
saveInputsTChem()
runTchem()
print('*-------------*')
print('Changing coverages ')
compareTchemAndCantera()<jupyter_output>*-------------*
Changing coverages 
-----------------------------
Gas phase
species name cantera  Tchem diff
-----------------------------
Ar 0.0 0.0 0.0 0.0
    
Ne 0.0 0.0 0.0 0.0
    
N2 0.0 0.0 0.0 0.0
    
CH4(2) 18289948.49921087 18290100.0 -151.50078912824392 -8.28328134082939e-06
    
O2(3) 315076.40871981625 315087.7 -11.291280183766503 -3.583664111713089e-05
    
CO2(4) 0.0 0.0 0.0 0.0
    
H2O(5) 497651.2875404093 497672.7 -21.41245959070511 -4.302703544993123e-05
    
H2(6) 9213972.47064621 9214129.0 -156.52935378998518 -1.6988259329909544e-05
    
CO(7) 9135958.043271543 9136126.0 -167.9567284565419 -1.8384139644800445e-05
    
C2H6(8) 3230637.2515739514 3230665.0 -27.748426048550755 -8.589149411631359e-06
    
CH2O(9) 128924.37016224995 128928.1 -3.7298377500555944 -2.893043220115509e-05
    
CH3(10) -18681821.433593906 -18681980.0 158.5664060935378 -8.487738021540101e-06
    
C3H8(11) 39900.46935646794 39901.35 -0.8806435320584569 -2.207100683931436e-05
    
H(12) -33153[...]<jupyter_text>Change gas mass fraction<jupyter_code>if make_random_sample:
    Ykk , Zkk = makeTChemRandomSample()
else:    
    Ykk , Zkk = getOriginalInput()
    
    Ykk[-1] *=2
    Ykk[-2] *=0.5
    Ykk[-3] *=1.5 

gas.Y  = Ykk[2:24]
gas.TP = Ykk[0], Ykk[1]
#surf.Y=Zkk[2:19]
surf.TP = Ykk[0], Ykk[1]
surf.coverages = Zkk
#surf.set_unnormalized_coverages =  cover
print('\nsurf.coverages after\n', surf.coverages)
saveInputsTChem()
runTchem()
print('*-------------*')
print('Changing molar gas fraction ')
compareTchemAndCantera()<jupyter_output>*-------------*
Changing molar gas fraction 
-----------------------------
Gas phase
species name cantera  Tchem diff
-----------------------------
Ar 0.0 0.0 0.0 0.0
    
Ne 0.0 0.0 0.0 0.0
    
N2 0.0 0.0 0.0 0.0
    
CH4(2) 45225078.51857975 45225180.0 -101.48142024874687 -2.2439191610702328e-06
    
O2(3) 2832471.7836432434 2832568.0 -96.21635675663128 -3.3969043332489547e-05
    
CO2(4) 0.0 0.0 0.0 0.0
    
H2O(5) 359369.2569348809 359389.5 -20.24306511908071 -5.63294292108822e-05
    
H2(6) 2062088.198043246 2062116.0 -27.801956753944978 -1.3482428530616088e-05
    
CO(7) 1564437.2004924177 1564457.0 -19.79950758232735 -1.265599384628243e-05
    
C2H6(8) 67712.22089399613 67715.55 -3.329106003875495 -4.916551192564234e-05
    
CH2O(9) -45489.63834380653 -45491.54 1.9016561934695346 -4.180416162241148e-05
    
CH3(10) -44153079.12832212 -44153130.0 50.87167788296938 -1.1521660298055542e-06
    
C3H8(11) 80.49122040573523 80.49255 -0.0013295942647602033 -1.65185005030122e-05
    
H[...]
{'fraction_non_alphanumeric': 0.1728892019724263, 'fraction_numerical': 0.1821475294354433, 'mean_word_length': 3.140833333333333, 'pattern_counts': {'":': 0, '<': 22, '<?xml version=': 0, '>': 24, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '6', 'dir': 'jupyter-structured-clean-dedup', 'id': '157967', 'license': 'BSD-2-Clause', 'n_tokens_mistral': 5081, 'n_tokens_neox': 4010, 'n_words': 623, 'path': 'src/example/SurfVerification/cantera/ComparingSurfCanteraTChem++.ipynb', 'repo_name': 'kyungjoo-kim/TChem', 'size': '9937'}