Skip to content

Commit

Permalink
Troe Tinkering
Browse files Browse the repository at this point in the history
  • Loading branch information
@tsikes
tsikes committed Mar 1, 2021
1 parent de2f6d6 commit 687737f
Show file tree
Hide file tree
Showing 3 changed files with 242 additions and 53 deletions.
254 changes: 208 additions & 46 deletions src/optimize/fit_coeffs.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -318,7 +318,7 @@ def nlopt_fit_fcn(x, grad):

return x

def fit_Troe(rates, T, P, X, rxnIdx, coefKeys, coefNames, mech, x0, bnds):
def fit_Troe_use_ct(rates, T, P, X, rxnIdx, coefKeys, coefNames, mech, x0, bnds):
def fit_rate_eqn(ln_k, P, X, mech, key, coefNames, rxnIdx):
rxn = mech.gas.reaction(rxnIdx)
def inner(temperatures, coeffs, scale_calc):
Expand Down Expand Up @@ -395,8 +395,8 @@ def inner(temperatures, coeffs, scale_calc):

return coeffs

def fit_Troe_no_ct(rates, T, M, x0=[], coefNames=default_Troe_coefNames, bnds=[]):
def fit_fcn_decorator(x0, alter_idx):
def fit_Troe(rates, T, M, x0=[], coefNames=default_Troe_coefNames, bnds=[], scipy_curvefit=True, Fit_LPL_HPL=False):
def fit_fcn_decorator(x0, M, alter_idx, s=[], jac=False):
def set_coeffs(*args):
coeffs = x0
for n, idx in enumerate(alter_idx):
Expand All @@ -410,7 +410,21 @@ def ln_Troe(T, *args):
k_inf = A_inf*T**n_inf*np.exp(-Ea_inf/(Ru*T))
P_r = k_0/k_inf*M
log_P_r = np.log10(P_r)
Fcent = (1-A)*np.exp(-T/T3)+A*np.exp(-T/T1)+np.exp(-T2/T)

if T3 == 0.0 or (T/T3 > max_ln_val).any():
exp_T_3 = 0
else:
exp_T_3 = np.exp(-T/T3)

if T1 == 0.0 or (T/T1 > max_ln_val).any():
exp_T_1 = 0
else:
exp_T_1 = np.exp(-T/T1)

Fcent = (1-A)*exp_T_3 + A*exp_T_1 + np.exp(-T2/T)
if (Fcent <= 0.0).any():
return np.ones_like(T)*np.inf

log_Fcent = np.log10(Fcent)
C = -0.4 - 0.67*log_Fcent
N = 0.75 - 1.27*log_Fcent
Expand All @@ -424,64 +438,189 @@ def ln_Troe(T, *args):

return ln_k

return ln_Troe
def ln_Troe_jac(T, *args): # TODO: currently not working, maybe use autodiff?
[Ea_0, ln_A_0, n_0, Ea_inf, ln_A_inf, n_inf, a, b, c, d, e] = set_coeffs(*args)
A_0, A_inf = np.exp(ln_A_0), np.exp(ln_A_inf)
k_0 = A_0*T**n_0*np.exp(-Ea_0/(Ru*T))
k_inf = A_inf*T**n_inf*np.exp(-Ea_inf/(Ru*T))
P_r = k_0/k_inf*M

if c == 0.0 or (-T/c > max_ln_val).any():
exp_neg_T_c = 0
else:
exp_neg_T_c = np.exp(-T/c)

abc_interior = a*np.exp(-b/T) + exp_neg_T_c
abc = 1/((1 + np.log10(P_r)**2)*abc_interior)

if (set([0, 1, 2, 3, 4, 5]) & set(alter_idx)): # if any arrhenius variable is being altered
ln_P_r_term = 1/(1 + P_r) - 2*np.log(abc_interior)*np.log10(P_r)/(1 + np.log10(P_r)**2)**2

jac = []
for n in alter_idx:
if n == 0: # dlnk_dEa_0
jac.append(-1/(Ru*T)*ln_P_r_term)
elif n == 1: # dlnk_dA_0
jac.append(1/A_0*ln_P_r_term)
elif n == 2: # dlnk_dn_0
jac.append(np.log(T)*ln_P_r_term)
elif n == 3: # dlnk_dEa_inf
jac.append(1/(Ru*T)*ln_P_r_term)
elif n == 4: # dlnk_dA_inf
jac.append(-1/A_inf*ln_P_r_term)
elif n == 5: # dlnk_dn_inf
jac.append(-np.log(T)*ln_P_r_term)
elif n == 6: # dlnk_da
jac.append(np.exp(-b/T)*abc)
elif n == 7: # dlnk_db
jac.append(-a/T*np.exp(-b/T)*abc)
elif n == 8: # dlnk_dc
if c == 0.0:
jac.append(np.zeros_like(T))
else:
jac.append(T/c**2*exp_neg_T_c*abc)
elif n == 9: # dlnk_d_d
jac.append(np.ones_like(T)/d)
elif n == 10:# dlnk_de
jac.append(np.log(T))

jac = np.vstack(jac).T
return jac

if not jac:
return ln_Troe
else:
return ln_Troe_jac

def nlopt_fit_fcn_decorator(fit_fcn, grad_fcn, x0, idx, T, ln_k_original):
def nlopt_fit_fcn(x, grad=[]):
x = x/s + x0[idx]

resid = fit_func(T, *x) - ln_k_original
loss = generalized_loss_fcn(resid).sum()

#s[:] = np.abs(np.sum(loss*grad_fcn(T, *x).T, axis=1))

#if len(grad) > 0:
# grad[:] = np.sum(loss*grad_fcn(T, *x).T, axis=1)

return loss
return nlopt_fit_fcn

ln_k = np.log(rates)

alter_idx = []
for n, coefName in enumerate(default_Troe_coefNames): # ['Ea_0', 'A_0', 'n_0', 'Ea_inf', 'A_inf', 'n_inf', 'A', 'T3', 'T1', 'T2']
alter_idx = {'low_rate': [], 'high_rate': [], 'falloff_parameters': [], 'all': []}
for n, coefName in enumerate(default_Troe_coefNames):
if coefName in coefNames:
alter_idx.append(n)

alter_idx['all'].append(n)
if coefName in ['Ea_0', 'A_0', 'n_0']:
alter_idx['low_rate'].append(n)
elif coefName in ['Ea_inf', 'A_inf', 'n_inf']:
alter_idx['high_rate'].append(n)
else:
alter_idx['falloff_parameters'].append(n)

if (set([0, 1, 2]) & set(alter_idx)) and len(x0) == 0:
a0 = np.polyfit(np.reciprocal(T[0:3]), ln_k[0:3], 1)
x0[0:3] = np.array([-a0[0]*Ru, np.exp(a0[1]), 0])
idx = alter_idx['low_rate']
a0 = np.polyfit(np.reciprocal(T[idx]), ln_k[idx], 1)
x0[idx] = np.array([-a0[0]*Ru, np.exp(a0[1]), 0])

if (set([3, 4, 5]) & set(alter_idx)) and len(x0) < 4:
a0 = np.polyfit(np.reciprocal(T[3:6]), ln_k[3:6], 1)
x0[3:6] = np.array([-a0[0]*Ru, np.exp(a0[1]), 0])
idx = alter_idx['high_rate']
a0 = np.polyfit(np.reciprocal(T[idx]), ln_k[idx], 1)
x0[idx] = np.array([-a0[0]*Ru, np.exp(a0[1]), 0])

# initial guesses for falloff
if len(x0) != 10:
x0 = [*x0[:6], 0.1, 100, 1000, 10000] # initial guesses for fitting Troe if none exist

x0 = np.array(x0)

if len(x0) < 7:
x0[6:9] = [0.1, 100, 1000, 10000] # initial guesses for fitting Troe if none exist
A_idx = [1, 4]
#A_idx = None
#if set(['A_0', 'A_inf']) & set(coefNames):
# A_idx = [i for i, coef in enumerate(coefNames) if coef in ['A_0', 'A_inf']]

x0[1] = np.log(x0[1])
x0[4] = np.log(x0[4])
# only valid initial guesses
bnds = deepcopy(bnds)
for n, val in enumerate(x0):
if val < bnds[0][n]:
x0[n] = bnds[0][n]
elif val > bnds[1][n]:
x0[n] = bnds[1][n]

x0 = np.array(x0)
for arrhenius_type in ['low_rate', 'high_rate']:
idx = alter_idx[arrhenius_type]
if len(idx) > 0:
x0[idx] = fit_arrhenius(rates[idx], T[idx], x0=x0[idx], bnds=[bnds[0][idx], bnds[1][idx]])

A_idx = None
if set(['A_0', 'A_inf']) & set(coefNames):
A_idx = np.argwhere(coefNames in ['A_0', 'A_inf'])
if A_idx is not None:
x0[A_idx] = np.log(x0[A_idx])
bnds[0][A_idx] = np.log(bnds[0][A_idx])
bnds[1][A_idx] = np.log(bnds[1][A_idx])

fit_func = fit_fcn_decorator(x0, alter_idx)
p0 = x0[alter_idx]
if not Fit_LPL_HPL:
idx = alter_idx['falloff_parameters']
T, M, ln_k = T[idx], M[idx], ln_k[idx]
p0 = x0[idx]

if len(bnds) > 0:
if A_idx is not None:
bnds[0][A_idx] = np.log(bnds[0][A_idx])
bnds[1][A_idx] = np.log(bnds[1][A_idx])
if scipy_curvefit:
fit_func = fit_fcn_decorator(x0, M, idx)
fit_func_jac = fit_fcn_decorator(x0, M, idx, jac=True)

# only valid initial guesses
for n, val in enumerate(p0):
if val < bnds[0][n]:
p0[n] = bnds[0][n]
elif val > bnds[1][n]:
p0[n] = bnds[1][n]
if len(bnds) == 0:
bnds = [np.ones_like(p0[idx]), np.ones_like(p0[idx])]*np.inf
else:
bnds = [bnds[0][idx], bnds[1][idx]]

with warnings.catch_warnings():
warnings.simplefilter('ignore', OptimizeWarning)
popt, _ = curve_fit(fit_func, T, ln_k, p0=p0, method='dogbox', bounds=bnds,
jac='2-point', x_scale='jac', max_nfev=len(p0)*1000)
else:
with warnings.catch_warnings():
warnings.simplefilter('ignore', OptimizeWarning)
popt, _ = curve_fit(fit_func, T, ln_k, p0=p0, method='dogbox',
jac='2-point', x_scale='jac', max_nfev=len(p0)*1000)
with warnings.catch_warnings():
warnings.simplefilter('ignore', OptimizeWarning)
#try:
x, _ = curve_fit(fit_func, T, ln_k, p0=p0, method='dogbox', bounds=bnds,
jac='3-point', x_scale='jac', max_nfev=len(p0)*1000)
#except:
# return

else:
s = np.ones_like(idx)
fit_func = fit_fcn_decorator(x0, M, idx, s=s)
fit_func_jac = lambda x: approx_fprime(x, lambda x: fit_func(T, x/s + x0[idx]), 1E-5)
nlopt_fit_fcn = nlopt_fit_fcn_decorator(fit_func, fit_func_jac, x0, idx, T, ln_k)


s[:] = fit_func_jac(np.zeros_like(p0))
print('s', s)
s[s == 0.0] = 1E-9
#s[s==0] = 10**(OoM(np.min(s[s!=0])) - 1) # TODO: MAKE THIS BETTER running into problem when s is zero, this is a janky workaround
s[:] = np.median(np.abs(s), axis=0)

opt = nlopt.opt(nlopt.LN_SBPLX, len(idx)) # nlopt.LN_SBPLX nlopt.LN_COBYLA nlopt.LD_MMA nlopt.LD_LBFGS

opt.set_min_objective(nlopt_fit_fcn)
#opt.set_maxeval(int(options['stop_criteria_val'])-1)
#opt.set_maxtime(options['stop_criteria_val']*60)

opt.set_xtol_rel(1E-2)
opt.set_ftol_rel(1E-2)
opt.set_lower_bounds((bnds[0][idx]-x0[idx])*s)
opt.set_upper_bounds((bnds[1][idx]-x0[idx])*s)

opt.set_initial_step(np.min(s[s != 1E-9]))
x = opt.optimize(np.zeros_like(p0))/s + x0[idx] # optimize!

x = np.array([*x0[:6], *x])

if A_idx is not None:
popt[A_idx] = np.exp(popt[A_idx])
print(f'x {x[-4:]}')
print(ln_k)
fit_k = fit_fcn_decorator(x0, M, idx)(T, *x[6:])
print(fit_k)
rss = np.sum((ln_k - fit_k)**2)
print(f'ln_k_resid {rss}')

return popt
if A_idx is not None:
x[A_idx] = np.exp(x[A_idx])

return x

def fit_generic(rates, T, P, X, rxnIdx, coefKeys, coefNames, mech, x0, bnds):
rxn = mech.gas.reaction(rxnIdx)
Expand All @@ -502,9 +641,28 @@ def fit_generic(rates, T, P, X, rxnIdx, coefKeys, coefNames, mech, x0, bnds):
coeffs[A_idx] = coeffs[A_idx]/mech.M(rxn)

elif type(rxn) is ct.FalloffReaction:
M = mech.M(rxn, [T, P, X])

falloff_coefNames = []
for key, coefName in zip(coefKeys, coefNames):
if coefName == 'activation_energy':
falloff_coefNames.append('Ea')
elif coefName == 'pre_exponential_factor':
falloff_coefNames.append('A')
elif coefName == 'temperature_exponent':
falloff_coefNames.append('n')

if key['coeffs'] == 'low_rate':
falloff_coefNames[-1] = f'{falloff_coefNames[-1]}_0'
elif key['coeffs'] == 'high_rate':
falloff_coefNames[-1] = f'{falloff_coefNames[-1]}_inf'

if rxn.falloff.type == 'Troe':
coeffs = fit_Troe(rates, T, P, X, rxnIdx, coefKeys, coefNames, mech, x0, bnds)
falloff_coefNames.extend(['A', 'T3', 'T1', 'T2'])
coeffs = fit_Troe(rates, T, M, x0=x0, coefNames=falloff_coefNames, bnds=bnds, scipy_curvefit=False)

elif rxn.falloff.type == 'SRI':
<<<<<<< Updated upstream
M = mech.M(rxn, [T, P, X])

SRI_coefNames = []
Expand All @@ -523,6 +681,10 @@ def fit_generic(rates, T, P, X, rxnIdx, coefKeys, coefNames, mech, x0, bnds):

SRI_coefNames.extend(['a', 'b', 'c', 'd', 'e'])
coeffs = fit_SRI(rates, T, M, x0, coefNames=SRI_coefNames, bnds=bnds)
=======
falloff_coefNames.extend(['a', 'b', 'c', 'd', 'e'])
coeffs = fit_SRI(rates, T, M, x0, coefNames=SRI_coefNames, bnds=bnds, scipy_curvefit=False)
>>>>>>> Stashed changes

return coeffs

Expand Down
34 changes: 28 additions & 6 deletions src/optimize/mech_optimize.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -14,7 +14,7 @@
from optimize.optimize_worker import Worker
from optimize.fit_fcn import initialize_parallel_worker, update_mech_coef_opt
from optimize.misc_fcns import rates, set_bnds
from optimize.fit_coeffs import fit_SRI, fit_Troe_no_ct
from optimize.fit_coeffs import fit_SRI, fit_Troe


default_arrhenius_coefNames = ['activation_energy', 'pre_exponential_factor', 'temperature_exponent']
Expand Down Expand Up @@ -302,6 +302,17 @@ def _set_rxn_rate_opt(self, rxn_coef_opt):

def _update_gas(self, rxn_coef_opt, rxn_rate_opt): # TODO: What happens if a second optimization is run?
mech = self.parent.mech
coef_opt = self.coef_opt

# delete any falloff coefficients with 5 indices
delete_idx = []
for i, idxDict in enumerate(coef_opt):
rxnIdx, coefName = idxDict['rxnIdx'], idxDict['coefName']
if coefName == 4:
delete_idx.append(i)

for i in delete_idx[::-1]:
del coef_opt[i]

rxns_changed = []
coeffs = []
Expand All @@ -320,6 +331,7 @@ def _update_gas(self, rxn_coef_opt, rxn_rate_opt): # TODO: What happens if a sec
if type(rxn) is ct.FalloffReaction:
lb = rxn_coef['coef_bnds']['lower']
ub = rxn_coef['coef_bnds']['upper']
<<<<<<< Updated upstream
if rxn.falloff.type == 'Troe':
rxn_coef['coef_x0'][6:] = fit_SRI(rates, T, M, x0=rxn_coef['coef_x0'],
coefNames=['a', 'b', 'c', 'd', 'e'], bnds=[lb, ub], scipy_curvefit=True)
Expand All @@ -328,19 +340,29 @@ def _update_gas(self, rxn_coef_opt, rxn_rate_opt): # TODO: What happens if a sec
coefNames=['a', 'b', 'c', 'd', 'e'], bnds=[lb, ub], scipy_curvefit=True)

mech.coeffs[rxnIdx]['falloff_type'] = 'SRI'
=======
if rxn.falloff.type == 'Troe': # This is for testing, it's not really needed
rxn_coef['coef_x0'] = fit_Troe(rates, T, M, x0=rxn_coef['coef_x0'], coefNames=['A', 'T3', 'T1', 'T2'],
bnds=[lb, ub], scipy_curvefit=True)
else:
rxn_coef['coef_x0'] = fit_Troe(rates, T, M, x0=rxn_coef['coef_x0'], coefNames=['A', 'T3', 'T1', 'T2'],
bnds=[lb, ub], scipy_curvefit=True)

mech.coeffs[rxnIdx]['falloff_type'] = 'Troe'
>>>>>>> Stashed changes
mech.coeffs[rxnIdx]['falloff_parameters'] = rxn_coef['coef_x0'][6:]

else:
lb = rxn_coef['coef_bnds']['lower']
ub = rxn_coef['coef_bnds']['upper']
rxn_coef['coef_x0'] = fit_SRI(rates, T, M, x0=rxn_coef['coef_x0'], bnds=[lb, ub], scipy_curvefit=True)
rxn_coef['coef_x0'] = fit_Troe(rates, T, M, x0=rxn_coef['coef_x0'], bnds=[lb, ub], scipy_curvefit=True)

rxn_coef['coefIdx'].extend(range(0, 5)) # extend to include falloff coefficients
rxn_coef['coefName'].extend(range(0, 5)) # extend to include falloff coefficients
rxn_coef['key'].extend([{'coeffs': 'falloff_parameters', 'coeffs_bnds': 'falloff_parameters'} for i in range(0, 5)])
rxn_coef['key'].extend([{'coeffs': 'falloff_parameters', 'coeffs_bnds': 'falloff_parameters'} for i in range(0, 4)])

# modify mech.coeffs from plog to falloff
mech.coeffs[rxnIdx] = {'falloff_type': 'SRI', 'high_rate': {}, 'low_rate': {}, 'falloff_parameters': rxn_coef['coef_x0'][6:],
mech.coeffs[rxnIdx] = {'falloff_type': 'Troe', 'high_rate': {}, 'low_rate': {}, 'falloff_parameters': rxn_coef['coef_x0'][6:],
'default_efficiency': 1.0, 'efficiencies': {}}

n = 0
Expand Down Expand Up @@ -378,8 +400,8 @@ def _update_gas(self, rxn_coef_opt, rxn_rate_opt): # TODO: What happens if a sec
continue

if type(rxn) is ct.FalloffReaction: # this means input reaction is falloff
if rxn.falloff.type == 'Troe':
rxn.falloff = ct.SriFalloff(mech.coeffs[rxnIdx]['falloff_parameters'])
if rxn.falloff.type == 'SRI':
rxn.falloff = ct.TroeFalloff(mech.coeffs[rxnIdx]['falloff_parameters'])

elif type(rxn) is ct.PlogReaction: # only need to generate a new mech if going from Plog -> SRI
#generate_new_mech = True
Expand Down
Loading

0 comments on commit 687737f

Please sign in to comment.