Skip to content

Commit

Permalink
More Progress
Browse files Browse the repository at this point in the history
  • Loading branch information
@tsikes
tsikes committed Apr 7, 2021
1 parent e5dec8f commit 36b155d
Show file tree
Hide file tree
Showing 3 changed files with 103 additions and 52 deletions.
138 changes: 91 additions & 47 deletions src/calculate/optimize/fit_coeffs.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -17,12 +17,9 @@
Ru = ct.gas_constant
# Ru = 1.98720425864083

min_pos_system_value = np.finfo(float).eps*(1E2)
max_pos_system_value = (np.finfo(float).max*(1E-20))**(1/3)
min_neg_system_value = -max_pos_system_value
min_ln_val = np.log(min_pos_system_value)
min_pos_system_value = (np.finfo(float).tiny*(1E20))**(1/2)
max_pos_system_value = (np.finfo(float).max*(1E-20))**(1/2)
max_ln_val = np.log(max_pos_system_value)
T_max = 6000

default_arrhenius_coefNames = ['activation_energy', 'pre_exponential_factor', 'temperature_exponent']
default_Troe_coefNames = ['Ea_0', 'A_0', 'n_0', 'Ea_inf', 'A_inf', 'n_inf', 'A', 'T3', 'T1', 'T2']
Expand All @@ -37,14 +34,24 @@
# (1, 1, 1)
]

troe_all_bnds = {'A': {'-': [min_neg_system_value, max_pos_system_value],
'+': [min_neg_system_value, max_pos_system_value]},
'T3': {'-': [min_neg_system_value, -T_max/max_ln_val],
'+': [-T_max/min_ln_val, max_pos_system_value]},
'T1': {'-': [min_neg_system_value, -T_max/max_ln_val],
'+': [-T_max/min_ln_val, max_pos_system_value]},
'T2': {'-': [-T_max*max_ln_val, -T_max*min_ln_val],
'+': [-T_max*max_ln_val, -T_max*min_ln_val]}}
#troe_min_pos_system_value = (np.finfo(float).tiny*(1E20))**(1/3)
#troe_max_pos_system_value = (np.finfo(float).max*(1E-20))**(1/3)
#troe_min_neg_system_value = -max_pos_system_value
#troe_min_ln_val = np.log(troe_min_pos_system_value)
#troe_max_ln_val = np.log(troe_max_pos_system_value)
#T_max = 6000
#troe_all_bnds = {'A': {'-': [troe_min_neg_system_value, troe_max_pos_system_value],
# '+': [troe_min_neg_system_value, troe_max_pos_system_value]},
# 'T3': {'-': [troe_min_neg_system_value, -T_max/troe_max_ln_val],
# '+': [-T_max/troe_min_ln_val, troe_max_pos_system_value]},
# 'T1': {'-': [troe_min_neg_system_value, -T_max/troe_max_ln_val],
# '+': [-T_max/troe_min_ln_val, troe_max_pos_system_value]},
# 'T2': {'-': [-T_max*troe_max_ln_val, -T_max*troe_min_ln_val],
# '+': [-T_max*troe_max_ln_val, -T_max*troe_min_ln_val]}}
troe_all_bnds = {'A': {'-': [-1E2, 1.0], '+': [-1E2, 1.0]},
'T3': {'-': [-1E30, -30], '+': [30.0, 1E30]},
'T1': {'-': [-1E30, -30], '+': [30.0, 1E30]},
'T2': {'-': [-1E30, 1E30], '+': [-1E30, 1E30]}}

def fit_arrhenius(rates, T, x0=[], coefNames=default_arrhenius_coefNames, bnds=[]):
def fit_fcn_decorator(x0, alter_idx, jac=False):
Expand Down Expand Up @@ -146,24 +153,24 @@ def ln_Troe(T, *x, grad_calc=True):
k_inf = A_inf*T**n_inf*np.exp(-Ea_inf/(Ru*T))
if ([k_0, k_inf] <= min_pos_system_value).any():
return np.ones_like(T)*max_pos_system_value
print(k_0)
print(k_inf)
P_r = k_0/k_inf*M
log_P_r = np.log10(P_r)

if T3 == 0.0 or (-T/T3 < -max_ln_val).any():
exp_T3 = 0
elif (-T/T3 > max_ln_val).any():
exp_T3 = max_pos_system_value
else:
exp_T3 = np.exp(-T/T3)

if T1 == 0.0 or (-T/T1 < -max_ln_val).any():
exp_T1 = 0
elif (-T/T1 > max_ln_val).any():
exp_T1 = max_pos_system_value
else:
exp_T1 = np.exp(-T/T1)

if (-T2/T > max_ln_val).any():
if T2 == 0.0:
exp_T2 = 0
elif (-T2/T > max_ln_val).any():
exp_T2 = max_pos_system_value
else:
exp_T2 = np.exp(-T2/T)
Expand All @@ -190,23 +197,23 @@ def ln_Troe_jac(T, *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))
if ([k_0, k_inf] <= min_pos_system_value).any():
return np.ones((len(T), len(alter_idx)))*max_pos_system_value
P_r = k_0/k_inf*M

if T3 == 0.0 or (-T/T3 < -max_ln_val).any():
exp_T3 = 0
elif (-T/T3 > max_ln_val).any():
exp_T3 = max_pos_system_value
else:
exp_T3 = np.exp(-T/T3)

if T1 == 0.0 or (-T/T1 < -max_ln_val).any():
exp_T1 = 0
elif (-T/T1 > max_ln_val).any():
exp_T1 = max_pos_system_value
else:
exp_T1 = np.exp(-T/T1)

if (-T2/T > max_ln_val).any():
if T2 == 0.0:
exp_T2 = 0
elif (-T2/T > max_ln_val).any():
exp_T2 = max_pos_system_value
else:
exp_T2 = np.exp(-T2/T)
Expand Down Expand Up @@ -305,7 +312,11 @@ def obj_fcn(x, grad=np.array([])):
#s[:] = np.abs(np.sum(loss*fit_func_jac(T, *x).T, axis=1))

if grad.size > 0:
grad[:] = np.sum(fit_func_jac(T, *x).T*resid, axis=1)
jac = fit_func_jac(T, *x)
if np.isfinite(jac).all():
grad[:] = np.sum(jac.T*resid, axis=1)
else:
grad[:] = np.ones_like(resid)*max_pos_system_value

return obj_val
return obj_fcn
Expand Down Expand Up @@ -338,10 +349,9 @@ def obj_fcn(x, grad=np.array([])):
x0 = [*x0[:6], 0.7, 200, 300, 400] # initial guesses for fitting Troe
x0 = np.array(x0)

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']]
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']]

# only valid initial guesses
bnds = deepcopy(bnds)
Expand Down Expand Up @@ -411,10 +421,14 @@ def obj_fcn(x, grad=np.array([])):

else:
#s[:] = calc_s(fit_func_jac(T, *p0))

opt = nlopt.opt(nlopt.GN_CRS2_LM, len(idx))
initial_step = calc_s(fit_func_jac(T, *p0))

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

opt.set_min_objective(obj_fcn)
#opt.set_maxeval(int(options['stop_criteria_val'])-1)
Expand All @@ -427,32 +441,66 @@ def obj_fcn(x, grad=np.array([])):
opt.set_lower_bounds(bnds[0][idx])
opt.set_upper_bounds(bnds[1][idx])

opt.set_initial_step(calc_s(fit_func_jac(T, *p0))*1E-3)
opt.set_initial_step(initial_step*1E-2)

local_opt = nlopt.opt(nlopt.LD_MMA, len(idx))
local_opt.set_xtol_rel(1E-3)
local_opt.set_ftol_rel(1E-3)
local_opt.set_initial_step(initial_step*1E-3)

opt.set_local_optimizer(local_opt)
x_fit = opt.optimize(p0) # optimize!

if HPL_LPL_defined: # Fit falloff
#x = np.array([*x0[:6], *(x_fit*s + x0[idx])])
x = np.array([*x0[:6], *x_fit])
else:
x = np.array(x_fit)
obj_val = obj_fcn(x)
if obj_val < HoF['obj_fcn']:
HoF['obj_fcn'] = obj_val

#obj_val = obj_fcn(x)
if opt.last_optimum_value() < HoF['obj_fcn']:
HoF['obj_fcn'] = opt.last_optimum_value()
HoF['coeffs'] = x

p0 = HoF['coeffs']
troe_bnds = []
for n, coef in enumerate(['A', 'T3', 'T1', 'T2']):
if HoF['coeffs'][-(4-n)] < 0:
troe_bnds.append(troe_all_bnds[coef]['-'])
else:
troe_bnds.append(troe_all_bnds[coef]['+'])

troe_bnds = np.array(troe_bnds)
bnds[:,-4:] = troe_bnds.T

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

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

opt.set_xtol_rel(1E-4)
opt.set_ftol_rel(1E-4)
#opt.set_lower_bounds((bnds[0][idx]-p0)*s)
#opt.set_upper_bounds((bnds[1][idx]-p0)*s)
opt.set_lower_bounds(bnds[0][idx])
opt.set_upper_bounds(bnds[1][idx])

opt.set_initial_step(calc_s(fit_func_jac(T, *p0))*1E-3)
x_fit = opt.optimize(p0) # optimize!

#if (x[-4:] == x0[-4:]).all():
# print('no change')
#else:
# print('fit values found')

#print(f'x {x[-4:]}')
print(ln_k)
fit_k = fit_func(T, *HoF['coeffs'], grad_calc=False)
print(fit_k)
#rss = np.sum((ln_k - fit_k)**2)
#print(ln_k)
#fit_k = fit_func(T, *HoF['coeffs'], grad_calc=False)
#print(fit_k)
#rss = 1/2*np.sum((ln_k - fit_k)**2)
#print(f'ln_k_resid {rss}')
print(obj_val)
print(opt.last_optimum_value())

if A_idx is not None:
HoF['coeffs'][A_idx] = np.exp(HoF['coeffs'][A_idx])
Expand Down Expand Up @@ -495,13 +543,9 @@ def fit_generic(rates, T, P, X, rxnIdx, coefKeys, coefNames, mech, x0, bnds):
elif key['coeffs'] == 'high_rate':
falloff_coefNames[-1] = f'{falloff_coefNames[-1]}_inf'

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

elif rxn.falloff.type == 'SRI':
falloff_coefNames.extend(['a', 'b', 'c', 'd', 'e'])
coeffs = fit_SRI(rates, T, M, x0, coefNames=SRI_coefNames, bnds=bnds, scipy_curvefit=True)
falloff_coefNames.extend(['A', 'T3', 'T1', 'T2'])
coeffs = fit_Troe(rates, T, M, x0=x0, coefNames=falloff_coefNames, bnds=bnds,
HPL_LPL_defined=True, scipy_curvefit=False)

return coeffs

Expand Down
2 changes: 1 addition & 1 deletion src/calculate/optimize/mech_optimize.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -342,7 +342,7 @@ def _update_gas(self, rxn_coef_opt, rxn_rate_opt): # TODO: What happens if a sec
ub = rxn_coef['coef_bnds']['upper']
if rxn.falloff.type == 'SRI':
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)
bnds=[lb, ub], scipy_curvefit=False)

mech.coeffs[rxnIdx]['falloff_type'] = 'Troe'
mech.coeffs[rxnIdx]['falloff_parameters'] = rxn_coef['coef_x0'][6:]
Expand Down
15 changes: 11 additions & 4 deletions src/calculate/optimize/misc_fcns.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,9 +7,11 @@


Ru = ct.gas_constant
min_pos_system_value = np.finfo(float).eps*(1E2)

min_pos_system_value = (np.finfo(float).tiny*(1E20))**(1/3)
max_pos_system_value = (np.finfo(float).max*(1E-20))**(1/3)
min_neg_system_value = -max_pos_system_value
T_min = 300
T_max = 6000

default_arrhenius_coefNames = ['activation_energy', 'pre_exponential_factor', 'temperature_exponent']
Expand Down Expand Up @@ -142,15 +144,20 @@ def set_bnds(mech, rxnIdx, keys, coefNames):
if coef_x0 > 0:
coef_bnds['lower'].append(0) # Ea shouldn't change sign
else:
coef_bnds['lower'].append(-Ru*T_max*np.log(max_pos_system_value))
coef_bnds['lower'].append(-1E-3*Ru*T_max*np.log(max_pos_system_value))
elif coefName == 'pre_exponential_factor':
coef_bnds['lower'].append(min_pos_system_value) # A should be positive
elif coefName == 'temperature_exponent':
coef_bnds['lower'].append(np.log(min_pos_system_value)/np.log(T_max))
elif not isinstance(coefName, int): # ints will be falloff, they will be taken care of below
coef_bnds['lower'].append(min_neg_system_value)

# set upper bnds
if coefName == 'activation_energy' and coef_x0 < 0: # Ea shouldn't change sign
coef_bnds['upper'].append(0)
if coefName == 'activation_energy':
if coef_x0 < 0: # Ea shouldn't change sign
coef_bnds['upper'].append(0)
else:
coef_bnds['upper'].append(-1E-3*Ru*T_min*np.log(min_pos_system_value))
elif coefName == 'temperature_exponent':
coef_bnds['upper'].append(np.log(max_pos_system_value)/np.log(T_max))
elif not isinstance(coefName, int):
Expand Down

0 comments on commit 36b155d

Please sign in to comment.