diff --git a/src/optimize/fit_coeffs.py b/src/optimize/fit_coeffs.py index dadf6bb..10b4d85 100644 --- a/src/optimize/fit_coeffs.py +++ b/src/optimize/fit_coeffs.py @@ -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): @@ -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): @@ -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 @@ -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) @@ -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 = [] @@ -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 diff --git a/src/optimize/mech_optimize.py b/src/optimize/mech_optimize.py index 32fb75a..09d3716 100644 --- a/src/optimize/mech_optimize.py +++ b/src/optimize/mech_optimize.py @@ -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'] @@ -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 = [] @@ -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) @@ -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 @@ -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 diff --git a/src/optimize/misc_fcns.py b/src/optimize/misc_fcns.py index e53f560..9879c47 100644 --- a/src/optimize/misc_fcns.py +++ b/src/optimize/misc_fcns.py @@ -161,8 +161,9 @@ def set_bnds(mech, rxnIdx, keys, coefNames): coef_bnds['exist'].append([True, True]) if type(rxn) in [ct.FalloffReaction, ct.PlogReaction]: - for SRI_coef in ['a', 'b', 'c', 'd', 'e']: + for SRI_coef in ['A', 'T3', 'T1', 'T2']: coef_bnds['exist'].append([False, False]) +<<<<<<< Updated upstream if SRI_coef == 'a': # this restriction isn't stricly necessary but can run into issues with log(-val) without coef_bnds['lower'].append(0.0) @@ -177,6 +178,10 @@ def set_bnds(mech, rxnIdx, keys, coefNames): coef_bnds['upper'].append(np.log(max_pos_system_value)) else: coef_bnds['upper'].append(max_pos_system_value) +======= + coef_bnds['lower'].append(min_neg_system_value) + coef_bnds['upper'].append(max_pos_system_value) +>>>>>>> Stashed changes coef_bnds['exist'] = np.array(coef_bnds['exist']) coef_bnds['lower'] = np.array(coef_bnds['lower'])