Skip to content

Commit

Permalink
Working
Browse files Browse the repository at this point in the history
  • Loading branch information
@tsikes
tsikes committed Apr 5, 2021
1 parent add6d9d commit 47bc21d
Show file tree
Hide file tree
Showing 2 changed files with 98 additions and 77 deletions.
153 changes: 84 additions & 69 deletions src/calculate/optimize/fit_coeffs.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -17,7 +17,8 @@
# Ru = 1.98720425864083

min_pos_system_value = np.finfo(float).eps*(1E2)
max_pos_system_value = np.finfo(float).max*(1E-20)
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)
max_ln_val = np.log(max_pos_system_value)

Expand Down Expand Up @@ -123,6 +124,8 @@ def ln_Troe(T, *x, grad_calc=True):
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_like(T)*max_pos_system_value
P_r = k_0/k_inf*M
log_P_r = np.log10(P_r)

Expand All @@ -140,11 +143,14 @@ def ln_Troe(T, *x, grad_calc=True):
else:
exp_T1 = np.exp(-T/T1)

exp_T2 = np.exp(-T2/T)
if (-T2/T > max_ln_val).any():
exp_T2 = max_pos_system_value
else:
exp_T2 = np.exp(-T2/T)

Fcent = (1-A)*exp_T3 + A*exp_T1 + exp_T2
if (Fcent <= 0.0).any():
return np.ones_like(T)*np.inf
return np.ones_like(T)*max_pos_system_value

log_Fcent = np.log10(Fcent)
C = -0.4 - 0.67*log_Fcent
Expand All @@ -153,7 +159,7 @@ def ln_Troe(T, *x, grad_calc=True):

e = np.exp(1)

ln_F = log_Fcent/np.log10(e)/(1+f1**2)
ln_F = log_Fcent/np.log10(e)/(1 + f1**2)

ln_k_calc = np.log(k_inf*P_r/(1 + P_r)) + ln_F

Expand All @@ -179,8 +185,11 @@ def ln_Troe_jac(T, *args):
exp_T1 = max_pos_system_value
else:
exp_T1 = np.exp(-T/T1)

exp_T2 = np.exp(-T2/T)

if (-T2/T > max_ln_val).any():
exp_T2 = max_pos_system_value
else:
exp_T2 = np.exp(-T2/T)

Fcent = (1-A)*exp_T3 + A*exp_T1 + exp_T2
if (Fcent <= 0.0).any():
Expand Down Expand Up @@ -210,14 +219,16 @@ def ln_Troe_jac(T, *args):
for n in alter_idx:
if n == 0: # dlnk_dEa_0
jac.append(1/(Ru*T)*(k_inf/Pr_denom + interior_term_2))
elif n == 1: # dlnk_dA_0
jac.append(1/A_0*(1-M_k0/Pr_denom - interior_term_2))
elif n == 1: # dlnk_dln_A_0
#jac.append(1/A_0*(1-M_k0/Pr_denom - interior_term_2)) # dlnk_dA_0
jac.append(1-M_k0/Pr_denom - interior_term_2)
elif n == 2: # dlnk_dn_0
jac.append(np.log(T)*(k_inf/Pr_denom - interior_term_2))
elif n == 3: # dlnk_dEa_inf
jac.append(1/(Ru*T)*(M_k0/Pr_denom - interior_term_2))
elif n == 4: # dlnk_dA_inf
jac.append(1/A_inf*(1-k_inf/Pr_denom + interior_term_2))
elif n == 4: # dlnk_dln_A_inf
#jac.append(1/A_inf*(1-k_inf/Pr_denom + interior_term_2)) # dlnk_dA_inf
jac.append(1-k_inf/Pr_denom + interior_term_2)
elif n == 5: # dlnk_dn_inf
jac.append(np.log(T)*(M_k0/Pr_denom + interior_term_2))
elif n == 6: # dlnk_dA
Expand Down Expand Up @@ -273,7 +284,7 @@ def obj_fcn(x, grad=[]):

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

print(obj_val)
return obj_val
return obj_fcn

Expand Down Expand Up @@ -318,87 +329,91 @@ def obj_fcn(x, grad=[]):
x0[n] = bnds[0][n]
elif val > bnds[1][n]:
x0[n] = bnds[1][n]

# Fit HPL and LPL (for falloff this is exact, otherwise a guess)
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]])

if HPL_LPL_defined:
# Fit HPL and LPL
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]])
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])

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 falloff
idx = alter_idx['falloff_parameters']
T, M, ln_k = T[idx], M[idx], ln_k[idx]
p0 = x0[idx]
#p0 = np.zeros_like(x0[idx])
s = 10**OoM(x0[idx])

if len(bnds) == 0:
bnds = [-np.ones_like(p0), np.ones_like(p0)]*np.inf
else:
bnds = [bnds[0][idx], bnds[1][idx]]
if HPL_LPL_defined: # Fit falloff
idx = alter_idx['falloff_parameters']
else:
idx = alter_idx['all']

fit_func = fit_fcn_decorator(ln_k, x0, M, idx, s=s)
fit_func_jac = fit_fcn_decorator(ln_k, x0, M, idx, s=s, jac=True)
T, M, ln_k = T[idx], M[idx], ln_k[idx]
p0 = x0[idx]
#p0 = np.zeros_like(x0[idx])
s = np.ones_like(x0[idx])

if scipy_curvefit:
# for testing scipy least_squares, curve_fit is a wrapper for this fcn
#obj_fcn = obj_fcn_decorator(fit_func, fit_func_jac, x0, idx, T, ln_k, return_sum=False)
#obj_fcn_jac = lambda x: fit_func_jac(T, *x)
if len(bnds) == 0:
bnds = [-np.ones_like(p0), np.ones_like(p0)]*np.inf
else:
bnds = [bnds[0][idx], bnds[1][idx]]

#res = least_squares(obj_fcn, p0, method='trf', bounds=bnds,
# jac=obj_fcn_jac, x_scale='jac', max_nfev=len(p0)*1000)
fit_func = fit_fcn_decorator(ln_k, x0, M, idx, s=s)
fit_func_jac = fit_fcn_decorator(ln_k, x0, M, idx, s=s, jac=True)

with warnings.catch_warnings():
warnings.simplefilter('ignore', OptimizeWarning)
try:
x_fit, _ = curve_fit(fit_func, T, ln_k, p0=p0, method='trf', bounds=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)
except:
return
if scipy_curvefit:
# for testing scipy least_squares, curve_fit is a wrapper for this fcn
#obj_fcn = obj_fcn_decorator(fit_func, fit_func_jac, x0, idx, T, ln_k, return_sum=False)
#obj_fcn_jac = lambda x: fit_func_jac(T, *x)

else:
s[:] = calc_s(fit_func_jac(T, *p0))
#res = least_squares(obj_fcn, p0, method='trf', bounds=bnds,
# jac=obj_fcn_jac, x_scale='jac', max_nfev=len(p0)*1000)

nlopt_fit_fcn = obj_fcn_decorator(fit_func, fit_func_jac, x0, idx, T, ln_k)
with warnings.catch_warnings():
warnings.simplefilter('ignore', OptimizeWarning)
try:
x_fit, _ = curve_fit(fit_func, T, ln_k, p0=p0, method='trf', bounds=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)
except:
return

else:
#s[:] = calc_s(fit_func_jac(T, *p0))
nlopt_fit_fcn = obj_fcn_decorator(fit_func, fit_func_jac, x0, idx, T, ln_k)

opt = nlopt.opt(nlopt.LN_SBPLX, len(idx)) # nlopt.LN_SBPLX nlopt.LN_COBYLA nlopt.LD_MMA nlopt.LD_LBFGS
opt = nlopt.opt(nlopt.GN_CRS2_LM, len(idx))
#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_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]-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_xtol_rel(1E-2)
opt.set_ftol_rel(1E-2)
#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])

print('p0 ', p0)
#opt.set_initial_step(np.min(s[s != 1E-9]))
x_fit = opt.optimize(p0) # optimize!
#opt.set_initial_step(calc_s(fit_func_jac(T, *p0))*1E-8)
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)

#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, *x_fit, grad_calc=False)
fit_k = fit_func(T, *x_fit, grad_calc=False)
#print(fit_k)
#rss = np.sum((ln_k - fit_k)**2)
#print(f'ln_k_resid {rss}')
rss = np.sum((ln_k - fit_k)**2)
print(f'ln_k_resid {rss}')

if A_idx is not None:
x[A_idx] = np.exp(x[A_idx])
Expand Down
22 changes: 14 additions & 8 deletions src/calculate/optimize/mech_optimize.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -14,7 +14,7 @@
from calculate.optimize.optimize_worker import Worker
from calculate.optimize.fit_fcn import update_mech_coef_opt
from calculate.optimize.misc_fcns import rates, set_bnds
from calculate.optimize.fit_coeffs import fit_SRI, fit_Troe
from calculate.optimize.fit_coeffs import fit_Troe


default_arrhenius_coefNames = ['activation_energy', 'pre_exponential_factor', 'temperature_exponent']
Expand Down Expand Up @@ -250,16 +250,22 @@ def _set_rxn_coef_opt(self, shock_conditions, min_T_range=1000, min_P_range=1E4)

rxn_coef['invT'].append(np.linspace(*invT_bnds, n_coef))
rxn_coef['T'].append(np.divide(10000, rxn_coef['invT'][-1]))
if coef_type == 'low_rate':
rxn_coef['P'].append(np.ones(n_coef)*1E-1) # Doesn't matter, will evaluate LPL and HPL
elif coef_type == 'high_rate':
rxn_coef['P'].append(np.ones(n_coef)*1E10) # Doesn't matter, will evaluate LPL and HPL
if type(rxn) is ct.PlogReaction:
if coef_type == 'low_rate':
rxn_coef['P'].append(np.ones(n_coef)*P_bnds[0])
elif coef_type == 'high_rate':
rxn_coef['P'].append(np.ones(n_coef)*P_bnds[1])
else: # is falloff type equation so these are dummy values
if coef_type == 'low_rate':
rxn_coef['P'].append(np.ones(n_coef)*1E-1) # will evaluate LPL
elif coef_type == 'high_rate':
rxn_coef['P'].append(np.ones(n_coef)*1E8) # will evaluate HPL

rxn_coef['T'].append(np.linspace(*T_bnds, 5))
rxn_coef['invT'].append(np.divide(10000, rxn_coef['T'][-1]))

if type(rxn) is ct.PlogReaction:
rxn_coef['P'].append(np.linspace(*P_bnds, 5))
rxn_coef['P'].append(np.linspace(*P_bnds, 7)[1:-1])
else:
rxn_coef['P'].append(np.ones(5)*P) # Doesn't matter, will evaluate median P for falloff

Expand Down Expand Up @@ -345,7 +351,7 @@ def _update_gas(self, rxn_coef_opt, rxn_rate_opt): # TODO: What happens if a sec
else:
lb = rxn_coef['coef_bnds']['lower']
ub = rxn_coef['coef_bnds']['upper']
rxn_coef['coef_x0'] = fit_Troe(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], HPL_LPL_defined=False, scipy_curvefit=False)

rxn_coef['coefIdx'].extend(range(0, 5)) # extend to include falloff coefficients
rxn_coef['coefName'].extend(range(0, 5)) # extend to include falloff coefficients
Expand Down Expand Up @@ -393,7 +399,7 @@ def _update_gas(self, rxn_coef_opt, rxn_rate_opt): # TODO: What happens if a sec
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
elif type(rxn) is ct.PlogReaction: # only need to generate a new mech if going from Plog -> Falloff
#generate_new_mech = True
for param_type in ['low_rate', 'high_rate', 'falloff_parameters']:
print(mech.coeffs[rxnIdx][param_type])
Expand Down

0 comments on commit 47bc21d

Please sign in to comment.