Skip to content

Commit

Permalink
Minor changes
Browse files Browse the repository at this point in the history
Consolidated bisymlog scaling factor
Added new % abs density gradient to options
  • Loading branch information
tsikes committed Nov 2, 2021
1 parent cf6328b commit 7bd4d64
Show file tree
Hide file tree
Showing 6 changed files with 63 additions and 38 deletions.
2 changes: 2 additions & 0 deletions src/calculate/convert_units.py
Original file line number Diff line number Diff line change
Expand Up @@ -98,6 +98,8 @@ def __init__(self, C=None, scaling_factor=2.0, base=10):
def set_C_heuristically(self, y, scaling_factor=None): # scaling factor: 1 looks loglike, 2 linear like
if scaling_factor is None:
scaling_factor = self.scaling_factor
else:
self.scaling_factor = scaling_factor

min_y = y.min()
max_y = y.max()
Expand Down
20 changes: 12 additions & 8 deletions src/calculate/optimize/fit_fcn.py
Original file line number Diff line number Diff line change
Expand Up @@ -74,8 +74,8 @@ def update_mech_coef_opt(mech, coef_opt, x):

def calculate_residuals(args_list):
def resid_func(t_offset, t_adjust, t_sim, obs_sim, t_exp, obs_exp, weights, obs_bounds=[],
loss_alpha=2, loss_c=1, loss_penalty=True, scale='Linear', DoF=1,
opt_type='Residual', verbose=False):
loss_alpha=2, loss_c=1, loss_penalty=True, scale='Linear',
bisymlog_scaling_factor=1.0, DoF=1, opt_type='Residual', verbose=False):

def calc_exp_bounds(t_sim, t_exp):
t_bounds = [max([t_sim[0], t_exp[0]])] # Largest initial time in SIM and Exp
Expand Down Expand Up @@ -116,7 +116,7 @@ def calc_exp_bounds(t_sim, t_exp):
obs_bounds = np.log10(np.abs(obs_bounds[ind])).squeeze()

elif scale == 'Bisymlog':
bisymlog = Bisymlog(C=None, scaling_factor=2.0)
bisymlog = Bisymlog(C=None, scaling_factor=bisymlog_scaling_factor)
bisymlog.set_C_heuristically(obs_exp)
obs_exp_bisymlog = bisymlog.transform(obs_exp)
obs_sim_interp_bisymlog = bisymlog.transform(obs_sim_interp)
Expand Down Expand Up @@ -204,7 +204,8 @@ def calc_density(x, data, dim=1):
# calculate time adjust with mse (loss_alpha = 2, loss_c =1)
time_adj_func = lambda t_adjust: resid_func(shock['opt_time_offset'], t_adjust*10**t_unc_OoM,
ind_var, obs_sim, obs_exp[:,0], obs_exp[:,1], weights, obs_bounds, scale=var['scale'],
DoF=len(coef_opt), opt_type=var['obj_fcn_type'])
bisymlog_scaling_factor= var['bisymlog_scaling_factor'], DoF=len(coef_opt),
opt_type=var['obj_fcn_type'])

res = minimize_scalar(time_adj_func, bounds=var['t_unc']/10**t_unc_OoM, method='bounded')
t_unc = res.x*10**t_unc_OoM
Expand All @@ -214,8 +215,9 @@ def calc_density(x, data, dim=1):
if loss_alpha == 3.0:
loss_alpha_fcn = lambda alpha: resid_func(shock['opt_time_offset'], t_unc,
ind_var, obs_sim, obs_exp[:,0], obs_exp[:,1], weights, obs_bounds,
loss_alpha=alpha, loss_c=var['loss_c'], loss_penalty=True,
scale=var['scale'], DoF=len(coef_opt), opt_type=var['obj_fcn_type'])
loss_alpha=alpha, loss_c=var['loss_c'], loss_penalty=True, scale=var['scale'],
bisymlog_scaling_factor= var['bisymlog_scaling_factor'], DoF=len(coef_opt),
opt_type=var['obj_fcn_type'])

res = minimize_scalar(loss_alpha_fcn, bounds=[-100, 2], method='bounded')
loss_alpha = res.x
Expand All @@ -227,8 +229,9 @@ def calc_density(x, data, dim=1):

output = resid_func(shock['opt_time_offset'], t_unc, ind_var, obs_sim, obs_exp[:,0], obs_exp[:,1],
weights, obs_bounds, loss_alpha=loss_alpha, loss_c=var['loss_c'],
loss_penalty=loss_penalty, scale=var['scale'], DoF=len(coef_opt),
opt_type=var['obj_fcn_type'], verbose=True)
loss_penalty=loss_penalty, scale=var['scale'],
bisymlog_scaling_factor= var['bisymlog_scaling_factor'],
DoF=len(coef_opt), opt_type=var['obj_fcn_type'], verbose=True)

output['shock'] = shock
output['independent_var'] = ind_var
Expand Down Expand Up @@ -262,6 +265,7 @@ def __init__(self, input_dict):
self.dist = self.parent.optimize.dist
self.opt_settings = {'obj_fcn_type': self.parent.optimization_settings.get('obj_fcn', 'type'),
'scale': self.parent.optimization_settings.get('obj_fcn', 'scale'),
'bisymlog_scaling_factor': self.parent.plot.signal.bisymlog.scaling_factor,
'loss_alpha': self.parent.optimization_settings.get('obj_fcn', 'alpha'),
'loss_c': self.parent.optimization_settings.get('obj_fcn', 'c'),
'bayes_dist_type': self.parent.optimization_settings.get('obj_fcn', 'bayes_dist_type'),
Expand Down
65 changes: 38 additions & 27 deletions src/calculate/reactors.py
Original file line number Diff line number Diff line change
Expand Up @@ -13,31 +13,32 @@


# list of all possible variables
all_var = {'Laboratory Time': {'SIM_name': 't_lab', 'sub_type': None},
'Shockwave Time': {'SIM_name': 't_shock', 'sub_type': None},
'Gas Velocity': {'SIM_name': 'vel', 'sub_type': None},
'Temperature': {'SIM_name': 'T', 'sub_type': None},
'Pressure': {'SIM_name': 'P', 'sub_type': None},
'Enthalpy': {'SIM_name': 'h', 'sub_type': ['total', 'species']},
'Entropy': {'SIM_name': 's', 'sub_type': ['total', 'species']},
'Density': {'SIM_name': 'rho', 'sub_type': None},
'Density Gradient': {'SIM_name': 'drhodz', 'sub_type': ['total', 'rxn']},
'% Density Gradient': {'SIM_name': 'perc_drhodz', 'sub_type': ['rxn']},
'Mole Fraction': {'SIM_name': 'X', 'sub_type': ['species']},
'Mass Fraction': {'SIM_name': 'Y', 'sub_type': ['species']},
'Concentration': {'SIM_name': 'conc', 'sub_type': ['species']},
'Net Production Rate': {'SIM_name': 'wdot', 'sub_type': ['species']},
'Creation Rate': {'SIM_name': 'wdotfor', 'sub_type': ['species']},
'Destruction Rate': {'SIM_name': 'wdotrev', 'sub_type': ['species']},
'Heat Release Rate': {'SIM_name': 'HRR', 'sub_type': ['total', 'rxn']},
'Delta Enthalpy (Heat of Reaction)':{'SIM_name': 'delta_h', 'sub_type': ['rxn']},
'Delta Entropy': {'SIM_name': 'delta_s', 'sub_type': ['rxn']},
'Equilibrium Constant': {'SIM_name': 'eq_con', 'sub_type': ['rxn']},
'Forward Rate Constant': {'SIM_name': 'rate_con', 'sub_type': ['rxn']},
'Reverse Rate Constant': {'SIM_name': 'rate_con_rev', 'sub_type': ['rxn']},
'Net Rate of Progress': {'SIM_name': 'net_ROP', 'sub_type': ['rxn']},
'Forward Rate of Progress': {'SIM_name': 'for_ROP', 'sub_type': ['rxn']},
'Reverse Rate of Progress': {'SIM_name': 'rev_ROP', 'sub_type': ['rxn']}}
all_var = {'Laboratory Time': {'SIM_name': 't_lab', 'sub_type': None},
'Shockwave Time': {'SIM_name': 't_shock', 'sub_type': None},
'Gas Velocity': {'SIM_name': 'vel', 'sub_type': None},
'Temperature': {'SIM_name': 'T', 'sub_type': None},
'Pressure': {'SIM_name': 'P', 'sub_type': None},
'Enthalpy': {'SIM_name': 'h', 'sub_type': ['total', 'species']},
'Entropy': {'SIM_name': 's', 'sub_type': ['total', 'species']},
'Density': {'SIM_name': 'rho', 'sub_type': None},
'Density Gradient': {'SIM_name': 'drhodz', 'sub_type': ['total', 'rxn']},
'% Density Gradient': {'SIM_name': 'perc_drhodz', 'sub_type': ['rxn']},
'\u00B1 % |Density Gradient|': {'SIM_name': 'perc_abs_drhodz', 'sub_type': ['rxn']},
'Mole Fraction': {'SIM_name': 'X', 'sub_type': ['species']},
'Mass Fraction': {'SIM_name': 'Y', 'sub_type': ['species']},
'Concentration': {'SIM_name': 'conc', 'sub_type': ['species']},
'Net Production Rate': {'SIM_name': 'wdot', 'sub_type': ['species']},
'Creation Rate': {'SIM_name': 'wdotfor', 'sub_type': ['species']},
'Destruction Rate': {'SIM_name': 'wdotrev', 'sub_type': ['species']},
'Heat Release Rate': {'SIM_name': 'HRR', 'sub_type': ['total', 'rxn']},
'Delta Enthalpy (Heat of Reaction)':{'SIM_name': 'delta_h', 'sub_type': ['rxn']},
'Delta Entropy': {'SIM_name': 'delta_s', 'sub_type': ['rxn']},
'Equilibrium Constant': {'SIM_name': 'eq_con', 'sub_type': ['rxn']},
'Forward Rate Constant': {'SIM_name': 'rate_con', 'sub_type': ['rxn']},
'Reverse Rate Constant': {'SIM_name': 'rate_con_rev', 'sub_type': ['rxn']},
'Net Rate of Progress': {'SIM_name': 'net_ROP', 'sub_type': ['rxn']},
'Forward Rate of Progress': {'SIM_name': 'for_ROP', 'sub_type': ['rxn']},
'Reverse Rate of Progress': {'SIM_name': 'rev_ROP', 'sub_type': ['rxn']}}

rev_all_var = {all_var[key]['SIM_name']:
{'name': key, 'sub_type': all_var[key]['sub_type']} for key in all_var.keys()}
Expand Down Expand Up @@ -74,16 +75,26 @@ def __call__(self, idx=None, units='CGS'): # units must be 'CGS' or 'SI'
parent = self.parent
if self.name == 'drhodz_tot':
self.value['SI'] = shock_fcns.drhodz(parent.states)

elif self.name == 'drhodz':
self.value['SI'] = shock_fcns.drhodz_per_rxn(parent.states)

elif self.name == 'perc_drhodz':
drhodz_tot = parent.drhodz_tot(units='SI')[:,None]
drhodz = parent.drhodz(units='SI').T
if not np.any(drhodz_tot):
self.value['SI'] = np.zeros_like(drhodz)
else:
self.value['SI'] = drhodz/np.abs(drhodz_tot)*100
#self.value['SI'] = drhodz/np.abs(drhodz).sum()

elif self.name == 'perc_abs_drhodz':
drhodz_tot = parent.drhodz_tot(units='SI')[:,None]
drhodz = parent.drhodz(units='SI').T
if not np.any(drhodz_tot):
self.value['SI'] = np.zeros_like(drhodz)
else:
self.value['SI'] = drhodz/np.abs(drhodz).sum(axis=1)[:,None]*100

else:
self.value['SI'] = getattr(parent.states, SIM_Dict[self.name])

Expand Down Expand Up @@ -285,7 +296,7 @@ def incident_shock_reactor(self, gas, details, t_end, **kwargs):
drhodz_tot=np.nan, drhodz=np.nan, perc_drhodz=np.nan)

reactor_vars = ['t_lab', 't_shock', 'z', 'A', 'vel', 'T', 'P', 'h_tot', 'h',
's_tot', 's', 'rho', 'drhodz_tot', 'drhodz', 'perc_drhodz',
's_tot', 's', 'rho', 'drhodz_tot', 'drhodz', 'perc_drhodz', 'perc_abs_drhodz',
'Y', 'X', 'conc', 'wdot', 'wdotfor', 'wdotrev',
'HRR_tot', 'HRR', 'delta_h', 'delta_s',
'eq_con', 'rate_con', 'rate_con_rev', 'net_ROP', 'for_ROP', 'rev_ROP']
Expand Down
9 changes: 7 additions & 2 deletions src/plot/base_plot.py
Original file line number Diff line number Diff line change
Expand Up @@ -29,6 +29,8 @@
from plot.custom_mpl_ticker_formatter import *
from timeit import default_timer as timer


bisymlog_scaling_factor = 1.5 # 1.0 log-like 2.0 linear-like
class Base_Plot(QtCore.QObject):
def __init__(self, parent, widget, mpl_layout):
super().__init__(parent)
Expand All @@ -40,6 +42,9 @@ def __init__(self, parent, widget, mpl_layout):
mpl.scale.register_scale(AbsoluteLogScale)
mpl.scale.register_scale(BiSymmetricLogScale)

# set layout to be tight
self.fig.tight_layout()

# Set plot variables
self.x_zoom_constraint = False
self.y_zoom_constraint = False
Expand All @@ -51,7 +56,7 @@ def __init__(self, parent, widget, mpl_layout):
self.autoScale = [True, True]

# Bisymlog
self.bisymlog = Bisymlog(C=None, scaling_factor=2.0)
self.bisymlog = Bisymlog(C=None, scaling_factor = bisymlog_scaling_factor)

# Connect Signals
self._draw_event_signal = self.canvas.mpl_connect('draw_event', self._draw_event)
Expand All @@ -71,7 +76,7 @@ def create_canvas(self):
scales = {'linear': True, 'log': 0, 'abslog': 0, 'bisymlog': 0}
for ax in self.ax:
ax.scale = {'x': scales, 'y': deepcopy(scales)}
ax.ticklabel_format(scilimits=(-4, 4), useMathText=True)
ax.ticklabel_format(scilimits=(-3, 4), useMathText=True)
ax.animateAxisLabels = False

# Get background
Expand Down
2 changes: 1 addition & 1 deletion src/plot/custom_mplscale.py
Original file line number Diff line number Diff line change
Expand Up @@ -124,7 +124,7 @@ def OoM(x):
max_OoM = raw_tick_OoM[-1]
min_dist = scale_ticklocs[2] - scale_ticklocs[1]

if vmin <= 0 <= vmax:
if vmin <= 0 and 0 <= vmax:
if min_dist > self.transform(10**zero_OoM):
min_dist = self.inverse_transform(min_dist)
zero_OoM = np.round(np.log10(np.abs(min_dist)))
Expand Down
3 changes: 3 additions & 0 deletions src/settings.py
Original file line number Diff line number Diff line change
Expand Up @@ -246,6 +246,9 @@ def shock_output(self):
parent.path[file] = parent.path['sim_dir'] / 'Sim {:d} - {:s}'.format(self.sim_num, file)

def sim_output(self, var_name): # takes variable name and creates path for it
if var_name == '\u00B1 % |Density Gradient|': # lots of invalid characters, replace
var_name = 'signed % Abs Density Gradient'

name = 'Sim {:d} - {:s}.txt'.format(self.sim_num, var_name)
self.parent.path[var_name] = self.parent.path['sim_dir'] / name

Expand Down

0 comments on commit 7bd4d64

Please sign in to comment.