Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Fix combinations of concentration ranges for RMG SA #131

Open
wants to merge 7 commits into
base: main
Choose a base branch
from
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
4 changes: 2 additions & 2 deletions t3/main.py
Original file line number Diff line number Diff line change
Expand Up @@ -805,8 +805,8 @@ def determine_reactions_based_on_sa(self) -> List[int]:
"""
reaction_keys, pdep_rxns_to_explore = list(), list()
if not os.path.isdir(self.paths['SA solver']):
self.logger.error("Could not find the path to the SA solver output folder.\n"
"Not performing refinement based on sensitivity analysis!")
self.logger.error(f"Could not find the path to the SA solver output folder:\n{self.paths['SA solver']}\n"
f"Not performing refinement based on sensitivity analysis!")
return reaction_keys

sa_dict_max = dict()
Expand Down
36 changes: 17 additions & 19 deletions t3/simulate/rmg_constantTP.py
Original file line number Diff line number Diff line change
Expand Up @@ -116,7 +116,6 @@ def set_up(self):
write_rmg_input_file(
rmg=self.generate_rmg_reactors_for_simulation(),
t3=self.t3,
iteration=1, # Does not matter for simulating or computing SA.
path=self.rmg_input_file,
walltime=self.t3['options']['max_RMG_walltime'],
)
Expand Down Expand Up @@ -236,17 +235,16 @@ def get_sa_coefficients(self) -> Optional[dict]:
sa_dict (Optional[dict]): An SA dictionary, structure is given in the docstring for T3/t3/main.py
"""
chem_to_rmg_rxn_index_map = get_chem_to_rmg_rxn_index_map(chem_annotated_path=self.paths['chem annotated'])
solver_path = os.path.join(self.paths['SA'], 'solver')
if not os.path.exists(solver_path):
self.logger.error("Could not find the path to RMG's SA solver output folder.")
if not os.path.exists(self.paths['SA solver']):
self.logger.error(f"Could not find the path to RMG's SA solver output folder:\n{self.paths['SA solver']}.")
return None
sa_files = list()
for file_ in os.listdir(solver_path):
for file_ in os.listdir(self.paths['SA solver']):
if 'sensitivity' in file_ and file_.endswith(".csv"):
sa_files.append(file_)
sa_dict = {'kinetics': dict(), 'thermo': dict(), 'time': list()}
for sa_file in sa_files:
df = pd.read_csv(os.path.join(solver_path, sa_file))
df = pd.read_csv(os.path.join(self.paths['SA solver'], sa_file))
for header in df.columns:
sa_type = None
if 'Time' in header:
Expand Down Expand Up @@ -320,7 +318,7 @@ def generate_rmg_reactors_for_simulation(self) -> dict:
v_vals = get_values_within_range(value_range=reactor['V'],
num=self.t3['options']['num_sa_per_volume_range'])
for t_val in t_vals:
for param in p_vals if ranged_p else v_vals:
for param in p_vals if (ranged_p or len(p_vals) == 1) else v_vals:
for species_list in species_lists:
new_reactor = {k: v for k, v in reactor.items() if k not in ['T', 'P', 'V']}
new_reactor['T'] = t_val
Expand All @@ -340,10 +338,10 @@ def get_species_concentration_lists_from_ranged_params(self) -> List[List[dict]]
species_lists = list()
spc_indices_w_ranges = [i for i, spc in enumerate(self.rmg['species'])
if isinstance(spc['concentration'], (list, tuple))]
species_list = [{'label': spc['label'], 'concentration': spc['concentration']} for spc in self.rmg['species']
if (isinstance(spc['concentration'], (float, int)) and spc['concentration'] > 0)
or spc['balance'] or not spc['reactive']]
species_vals = [get_values_within_range(value_range=self.rmg['species'][spc_indices_w_ranges.index(species_index)]['concentration'],
species_list_wo_ranges = [{'label': spc['label'], 'concentration': spc['concentration']} for spc in self.rmg['species']
if (isinstance(spc['concentration'], (float, int)) and spc['concentration'] > 0)
or (spc['balance'] and spc['concentration'] == 0)]
species_vals = [get_values_within_range(value_range=self.rmg['species'][species_index]['concentration'],
num=self.t3['options']['num_sa_per_concentration_range'])
for species_index in spc_indices_w_ranges]

Expand All @@ -353,13 +351,13 @@ def get_species_concentration_lists_from_ranged_params(self) -> List[List[dict]]
if spc['concentration'] > 0 or spc['balance'] or not spc['reactive']]]

# 2. Only two ranged concentrations and modify_concentration_ranges_in_reverse is True
if len(spc_indices_w_ranges) == 2 and self.t3['options']['modify_concentration_ranges_in_reverse']:
elif len(spc_indices_w_ranges) == 2 and self.t3['options']['modify_concentration_ranges_in_reverse']:
spc_0_vals = get_values_within_range(value_range=self.rmg['species'][spc_indices_w_ranges[0]]['concentration'],
num=self.t3['options']['num_sa_per_concentration_range'])
spc_1_vals = get_values_within_range(value_range=self.rmg['species'][spc_indices_w_ranges[1]]['concentration'],
num=self.t3['options']['num_sa_per_concentration_range'])
for val_0, val_1 in zip(spc_0_vals, spc_1_vals[::-1]):
new_species_list = species_list
new_species_list = species_list_wo_ranges
new_species_list.append({'label': self.rmg['species'][spc_indices_w_ranges[0]]['label'],
'concentration': val_0})
new_species_list.append({'label': self.rmg['species'][spc_indices_w_ranges[1]]['label'],
Expand All @@ -368,17 +366,17 @@ def get_species_concentration_lists_from_ranged_params(self) -> List[List[dict]]

# 3. No combinations, modify_concentration_ranges_together is True
elif self.t3['options']['modify_concentration_ranges_together']:
for point_number in range(self.t3['options']['num_sa_per_concentration_range']):
new_species_list = species_list
for i, spc_index in enumerate(spc_indices_w_ranges):
new_species_list.append({'label': self.rmg['species'][spc_index]['label'],
'concentration': species_vals[i][point_number]})
for point_number in range(self.t3['options']['num_sa_per_concentration_range']):
new_species_list = [entry for entry in species_list_wo_ranges]
for i, spc_index in enumerate(spc_indices_w_ranges):
new_species_list.append({'label': self.rmg['species'][spc_index]['label'],
'concentration': species_vals[i][point_number]})
species_lists.append(new_species_list)

# 4. Combinations (products)
else:
for vals in itertools.product(*species_vals):
new_species_list = species_list
new_species_list = species_list_wo_ranges
for i, val in enumerate(vals):
new_species_list.append({'label': self.rmg['species'][spc_indices_w_ranges[i]]['label'],
'concentration': val})
Expand Down
22 changes: 18 additions & 4 deletions t3/utils/writer.py
Original file line number Diff line number Diff line change
Expand Up @@ -18,8 +18,9 @@

def write_rmg_input_file(rmg: dict,
t3: dict,
iteration: int,
path: str,
iteration: int = 1,
run_sa: bool = False,
walltime: str = '00:00:00:00',
):
"""
Expand All @@ -29,8 +30,10 @@ def write_rmg_input_file(rmg: dict,
Args:
rmg (dict): The arguments to write in a keyword argument dictionary format.
t3 (dict): The T3 arguments in a keyword argument dictionary format. Includes atol and rtol for SA.
iteration (int): The T3 iteration, used to determine ``core_tolerance`` and ``tolerance_interrupt_simulation``.
path (str): The path where the RMG input file should be saved.
iteration (int, optional): The T3 iteration, used to determine ``core_tolerance`` and
``tolerance_interrupt_simulation``. Does not matter for simulating or computing SA.
run_sa (bool, optional): Whether to add sensitivity analysis information into the RMG input file.
walltime (str, optional): The time cap for an RMG run. Should pass here t3['options']['max_RMG_walltime']
"""
rmg = rmg.copy()
Expand Down Expand Up @@ -100,7 +103,7 @@ def write_rmg_input_file(rmg: dict,
pressure=${pressure},
initialMoleFractions={${concentrations()} },
${termination}
nSims=${conditions_per_iteration},${balance}${constant}
nSims=${conditions_per_iteration},${balance}${constant}${sensitivity}
)
<%def name="concentrations()">
% for spc in species_list:
Expand All @@ -118,7 +121,7 @@ def write_rmg_input_file(rmg: dict,
temperature=${temperature},
initialConcentrations={${concentrations()} },
${termination}
nSims=${conditions_per_iteration},${constant}
nSims=${conditions_per_iteration},${constant}${sensitivity}
)
<%def name="concentrations()">
% for spc in species_list:
Expand Down Expand Up @@ -166,6 +169,16 @@ def write_rmg_input_file(rmg: dict,
constant = '\n constantSpecies=['
constant += f"'{spc['label']}', "
constant += '],' if constant else ''
sensitivity = ''
observables = [spc['label'] for spc in species if spc['observable']]
if run_sa and len(observables):
sensitivity = '\n sensitivity=['
for i, observable in enumerate(observables):
sensitivity += f"'{observable}'"
if i < len(observables) - 1:
sensitivity += ', '
sensitivity += f"],\n sensitivityThreshold={t3['sensitivity']['SA_threshold']}"


if reactor['type'] == 'gas batch constant T P':
if isinstance(reactor['P'], float):
Expand All @@ -188,6 +201,7 @@ def write_rmg_input_file(rmg: dict,
conditions_per_iteration=reactor['conditions_per_iteration'],
balance=balance,
constant=constant,
sensitivity=sensitivity,
)

elif reactor['type'] == 'liquid batch constant T V':
Expand Down