Skip to content

Commit

Permalink
Semi working
Browse files Browse the repository at this point in the history 
Troe is working, but need to implement constraints on Fcent fitting.

Should be working for Plog -> Troe after implementing those constraints
  • Loading branch information
@tsikes
tsikes committed Jun 24, 2021
1 parent afdcdbc commit d930224
Show file tree
Hide file tree
Showing 2 changed files with 77 additions and 38 deletions.
95 changes: 61 additions & 34 deletions src/calculate/optimize/fit_coeffs.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -576,33 +576,42 @@ def ln_Troe(M, *x):
return ln_k_calc
return ln_Troe

def fit_Fcent_decorator(Fcent, jac=False):
def Fcent_calc(T, *x):
def fit_ln_Fcent_decorator(Fcent, jac=False):
def set_coeffs(x):
A = x[0]
[T3, T1, T2] = np.exp(x[1:])
T3, T1 = 1000/np.exp(x[1]), 1000/np.exp(x[2])
T2 = 100*np.exp(x[3])

res = (1-A)*np.exp(-T/T3) + A*np.exp(-T/T1) + np.exp(-T2/T)
#A = np.log(x[0])
#T3, T1 = 1000/x[1], 1000/x[2]
#T2 = x[3]*100

return res
return [A, T3, T1, T2]

def Fcent_calc_jac(T, *x):
A = x[0]
[T3, T1, T2] = np.exp(x[1:])
def ln_Fcent_calc(T, *x):
[A, T3, T1, T2] = set_coeffs(x)

Fcent_fit = (1-A)*np.exp(-T/T3) + A*np.exp(-T/T1) + np.exp(-T2/T)

return np.log(Fcent_fit)

def ln_Fcent_calc_jac(T, *x): # Not for ln(Fcent) need to redo if using
[A, T3, T1, T2] = set_coeffs(x)

jac = []
jac.append(np.exp(-T/T1) - np.exp(-T/T3)) # dFcent/dA
jac.append((1-A)*T/(T3**2)*np.exp(-T/T3)) # dFcent/dT3
jac.append(A*T/(T1**2)*np.exp(-T/T1)) # dFcent/dT1
jac.append(-1/T*np.exp(-T2/T)) # dFcent/dT2
jac.append(1/A*(np.exp(-T/T1) - np.exp(-T/T3))) # dln_Fcent/dA
jac.append((1-A)*T/(T3**3)*np.exp(-T/T3)) # dln_Fcent/dT3
jac.append(A*T/(T1**3)*np.exp(-T/T1)) # dln_Fcent/dT1
jac.append(-1/(T2*T)*np.exp(-T2/T)) # dln_Fcent/dT2

jac = np.vstack(jac).T

return jac

if not jac:
return Fcent_calc
return ln_Fcent_calc
elif jac:
return Fcent_calc_jac
return ln_Fcent_calc_jac

def calc_s(grad):
x = np.abs(grad)
Expand Down Expand Up @@ -673,9 +682,9 @@ def obj_fcn(x, grad=np.array([])):
if len(res) == 0:
p0 = [27, 14, -0.5] # ln(k_0), ln(k_inf), ln(Fcent)
else:
p0 = np.log10(res[-1][1:])
p0 = np.log(res[-1][1:])

p_bnds = [[-27.631, -27.631, -18.421], [69.078, 69.078, 0]] # ln(k_0), ln(k_inf), ln(Fcent)
p_bnds = np.log([[1E-12, 1E-12, 1E-8], [1E30, 1E30, 1]]) # ln(k_0), ln(k_inf), ln(Fcent)

x_fit, _ = curve_fit(fit_func, M[idx], ln_k[idx], p0=p0, method='trf', bounds=p_bnds, # dogbox
#jac=fit_func_jac, x_scale='jac', max_nfev=len(p0)*1000)
Expand All @@ -697,27 +706,33 @@ def obj_fcn(x, grad=np.array([])):
x[idx] = fit_arrhenius(res[:, res_idx+1], res[:,0], bnds=[bnds[0][idx], bnds[1][idx]], loss='huber') # 'soft_l1', 'huber', 'cauchy', 'arctan'

T = res[:,0]
cmp = np.array([10000/T, np.ln(res[:, res_idx+1]), np.ln(x[idx[1]]*T**x[idx[2]]*np.exp(-x[idx[0]]/Ru/T))]).T
cmp = np.array([10000/T, np.log(res[:, res_idx+1]), np.log(x[idx[1]]*T**x[idx[2]]*np.exp(-x[idx[0]]/Ru/T))]).T
for entry in cmp:
print(*entry)
print('')

# Fit falloff
idx = alter_idx['falloff_parameters']
fit_func = fit_Fcent_decorator(res[:,3])
fit_func_jac = fit_Fcent_decorator(res[:,3], jac=True)
p0 = [0.6, 5.3, 6.2, 6.7]
p_bnds = [[-100, 3.4, 3.4, 3.4], [1, 69.078, 69.078, 69.078]]
fit_func = fit_ln_Fcent_decorator(res[:,3])
fit_func_jac = fit_ln_Fcent_decorator(res[:,3], jac=True)
p0 = [0.6, np.log(1000/200), np.log(1000/600), np.log(1200/100)]
p_bnds = [[-10, np.log(1E-37), np.log(1E-37), np.log(30/100)],
[1, np.log(1000/30), np.log(1000/30), np.log(1E38)]]

#p0 = [np.exp(0.6), 1000/200, 1000/600, 1200/100]
#p_bnds = [[np.exp(-100), 1E-37, 1E-37, 30/100], [np.exp(1), 1000/30, 1000/30, 1E38]]

x_fit, _ = curve_fit(fit_func, res[:,0], res[:,3], p0=p0, method='trf', bounds=p_bnds, # dogbox
x_fit, _ = curve_fit(fit_func, res[:,0], np.log(res[:,3]), p0=p0, method='trf', bounds=p_bnds, # dogbox
#jac=fit_func_jac, x_scale='jac', max_nfev=len(p0)*1000)
jac='2-point', x_scale='jac', max_nfev=len(p0)*1000, loss='huber')

cmp = np.array([res[:,0], res[:,3], fit_func(res[:,0], *x_fit)]).T
cmp = np.array([res[:,0], res[:,3], np.exp(fit_func(res[:,0], *x_fit))]).T
for entry in cmp:
print(*entry)

x[idx] = [x_fit[0], *np.exp(x_fit[1:])]
x[idx] = [x_fit[0], 1000/np.exp(x_fit[1]), 1000/np.exp(x_fit[2]), 100*np.exp(x_fit[3])]

#x[idx] = [np.log(x_fit[0]), 1000/x_fit[1], 1000/x_fit[2], 100*x_fit[3]]

else:
# only valid initial guesses
Expand All @@ -744,11 +759,11 @@ def obj_fcn(x, grad=np.array([])):
#fit_func = lambda x: (fit_const_T_decorator(ln_k[idx], T[idx])(M[idx], [ln_k_0[idx], ln_k_inf[idx], x[0]]) - ln_k[idx])**2 # only fitting Fcent
fit_func = lambda M, x: fit_const_T_decorator(ln_k[idx], T[idx])(M, [ln_k_0[idx], ln_k_inf[idx], x])
if len(res) == 0:
p0 = [-0.5] # log(k_0), log(k_inf), log(Fcent)
p0 = [-0.5] # ln(k_0), ln(k_inf), ln(Fcent)
else:
p0 = np.log(res[-1][1:])

p_bnds = [[-18.421], [0]] # ln(k_0), ln(k_inf), ln(Fcent) not sure if lower Fcent bnd should be -2 or -1
p_bnds = [[np.log(1E-8)], [0]] # ln(Fcent)

with warnings.catch_warnings():
warnings.simplefilter('ignore', OptimizeWarning)
Expand All @@ -760,29 +775,41 @@ def obj_fcn(x, grad=np.array([])):

res.append([T[idx], *np.exp(x_fit)])

#cmp = np.array([T[idx], M[idx], ln_k[idx], fit_func(M[idx], *x_fit)])
#print(*cmp)

res = np.array(res)

# Fit falloff
idx = alter_idx['falloff_parameters']
fit_func = fit_Fcent_decorator(res[:,1])
fit_func_jac = fit_Fcent_decorator(res[:,1], jac=True)
p0 = [0.6, 5.3, 6.2, 6.7]
p_bnds = [[-100, 3.4, 3.4, 3.4], [1, 69.078, 69.078, 69.078]]
fit_func = fit_ln_Fcent_decorator(res[:,1])
fit_func_jac = fit_ln_Fcent_decorator(res[:,1], jac=True)
p0 = [0.6, np.log(1000/200), np.log(1000/600), np.log(1200/100)]
p_bnds = [[-10, np.log(1E-37), np.log(1E-37), np.log(30/100)],
[1, np.log(1000/30), np.log(1000/30), np.log(1E38)]]

#p0 = [np.exp(0.6), 1000/200, 1000/600, 1200/100]
#p_bnds = [[np.exp(-100), 1E-37, 1E-37, 30/100], [np.exp(1), 1000/30, 1000/30, 1E38]]

with warnings.catch_warnings():
warnings.simplefilter('ignore', OptimizeWarning)
x_fit, _ = curve_fit(fit_func, res[:,0], res[:,1], p0=p0, method='trf', bounds=p_bnds, # dogbox
x_fit, _ = curve_fit(fit_func, res[:,0], np.log(res[:,1]), p0=p0, method='trf', bounds=p_bnds, # dogbox
#jac=fit_func_jac, x_scale='jac', max_nfev=len(p0)*1000)
jac='2-point', x_scale='jac', max_nfev=len(p0)*1000, loss='huber')

#cmp = np.array([res[:,0], res[:,1], fit_func(res[:,0], *x_fit)]).T
#for entry in cmp:
# print(*entry)
#print('')

x[idx] = [x_fit[0], 1000/np.exp(x_fit[1]), 1000/np.exp(x_fit[2]), 100*np.exp(x_fit[3])]

#x[idx] = [np.log(x_fit[0]), 1000/x_fit[1], 1000/x_fit[2], 100*x_fit[3]]

x[idx] = [x_fit[0], *np.exp(x_fit[1:])]
#Fcent = (1-A)*np.exp(-T/T3) + A*np.exp(-T/T1) + np.exp(-T2/T)
#fit_func = lambda M, x: fit_const_T_decorator(ln_k[idx], T[idx])(M, [ln_k_0[idx], ln_k_inf[idx], Fcent])

Fcent = (1-A)*np.exp(-T/T3) + A*np.exp(-T/T1) + np.exp(-T2/T)
fit_func = lambda M, x: fit_const_T_decorator(ln_k[idx], T[idx])(M, [ln_k_0[idx], ln_k_inf[idx], Fcent])
print(x[-4:])

return x

Expand Down
20 changes: 16 additions & 4 deletions src/calculate/optimize/mech_optimize.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -16,7 +16,7 @@
from calculate.optimize.misc_fcns import rates, set_bnds
from calculate.optimize.fit_coeffs import fit_Troe


Ru = ct.gas_constant
default_arrhenius_coefNames = ['activation_energy', 'pre_exponential_factor', 'temperature_exponent']

class Multithread_Optimize:
Expand Down Expand Up @@ -295,12 +295,24 @@ def _set_rxn_rate_opt(self):
i = 0
for rxn_coef in self.rxn_coef_opt: # LB and UB are off for troe arrhenius parts
rxnIdx = rxn_coef['rxnIdx']
rxn = mech.gas.reaction(rxnIdx)
rate_bnds_val = mech.rate_bnds[rxnIdx]['value']
rate_bnds_type = mech.rate_bnds[rxnIdx]['type']
for T, P in zip(rxn_coef['T'], rxn_coef['P']):
bnds = mech.rate_bnds[rxnIdx]['limits'](np.exp(rxn_rate_opt['x0'][i]))
for n, (T, P) in enumerate(zip(rxn_coef['T'], rxn_coef['P'])):
if type(rxn) is ct.FalloffReaction: # if falloff, change arrhenius rates to LPL/HPL
coefName = rxn_coef['coefName'][n]
if coefName in default_arrhenius_coefNames:
key = rxn_coef['key'][n]
x = []
for ArrheniusCoefName in default_arrhenius_coefNames:
x.append(mech.coeffs_bnds[rxnIdx][key['coeffs_bnds']][ArrheniusCoefName]['resetVal'])

rxn_rate_opt['x0'][i] = np.log(x[1]) + x[2]*np.log(T) - x[0]/(Ru*T)

ln_rate = rxn_rate_opt['x0'][i]
bnds = mech.rate_bnds[rxnIdx]['limits'](np.exp(ln_rate))
bnds = np.sort(np.log(bnds)) # operate on ln and scale
scaled_bnds = np.sort(bnds/rxn_rate_opt['x0'][i])
scaled_bnds = np.sort(bnds/ln_rate)
lb.append(scaled_bnds[0])
ub.append(scaled_bnds[1])

Expand Down

0 comments on commit d930224

Please sign in to comment.