Skip to content

Commit a6b61ca

Browse files
alongdalongd
committed
Auto stash before rebase of "idt"
1 parent 405944f commit a6b61ca

7 files changed

Lines changed: 491 additions & 16 deletions

File tree

‎t3/main.py‎

Lines changed: 104 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -158,7 +158,7 @@ def __init__(self,
158158
clean_dir: bool = False,
159159
):
160160

161-
self.sa_dict = None
161+
self.sa_dict, self.sa_dict_idt = None, None
162162
self.sa_observables = list()
163163
self.t0 = datetime.datetime.now() # initialize the timer as datetime object
164164

@@ -313,9 +313,12 @@ def execute(self):
313313
sa_rtol=self.t3['sensitivity']['rtol'],
314314
)
315315
simulate_adapter.simulate()
316-
self.sa_dict = simulate_adapter.get_sa_coefficients()
317-
save_yaml_file(path=self.paths['SA dict'], content=self.sa_dict)
318-
if self.t3['sensitivity']['global_observables'] == 'IDT':
316+
self.sa_dict = simulate_adapter.get_sa_coefficients(
317+
top_SA_species=self.t3['sensitivity']['top_SA_species'],
318+
top_SA_reactions=self.t3['sensitivity']['top_SA_reactions'],
319+
max_workers=self.t3['sensitivity']['max_workers'],
320+
save_yaml=True)
321+
if self.t3['sensitivity']['global_observables'].lower() == 'idt':
319322
simulate_adapter = simulate_factory(simulate_method='CanteraIDT',
320323
t3=self.t3,
321324
rmg=self.rmg,
@@ -327,8 +330,11 @@ def execute(self):
327330
sa_rtol=self.t3['sensitivity']['rtol'],
328331
)
329332
simulate_adapter.simulate()
330-
self.sa_dict_idt = simulate_adapter.get_sa_coefficients()
331-
save_yaml_file(path=self.paths['SA idt dict'], content=self.sa_dict_idt)
333+
self.sa_dict_idt = simulate_adapter.get_sa_coefficients(
334+
top_SA_species=self.t3['sensitivity']['top_SA_species'],
335+
top_SA_reactions=self.t3['sensitivity']['top_SA_reactions'],
336+
max_workers=self.t3['sensitivity']['max_workers'],
337+
save_yaml=True)
332338

333339
additional_calcs_required = self.determine_species_and_reactions_to_calculate()
334340

@@ -395,6 +401,7 @@ def set_paths(self,
395401
'SA input': os.path.join(iteration_path, 'SA', 'input.py'),
396402
'SA dict': os.path.join(iteration_path, 'SA', 'sa.yaml'),
397403
'SA IDT dict': os.path.join(iteration_path, 'SA', 'sa_idt.yaml'),
404+
'SA IDT dict top X': os.path.join(iteration_path, 'SA', 'sa_idt_top_x.yaml'),
398405
'PDep SA': os.path.join(iteration_path, 'PDep_SA'),
399406
'ARC': os.path.join(iteration_path, 'ARC'),
400407
'ARC input': os.path.join(iteration_path, 'ARC', 'input.yml'),
@@ -691,6 +698,7 @@ def determine_species_and_reactions_to_calculate(self) -> bool:
691698
bool: Whether additional calculations are required.
692699
"""
693700
species_keys, reaction_keys, coll_vio_spc_keys, coll_vio_rxn_keys = list(), list(), list(), list()
701+
rxn_idt_keys = None
694702

695703
self.rmg_species, self.rmg_reactions = self.load_species_and_reactions_from_chemkin_file()
696704
self.logger.info(f'This RMG model has {len(self.rmg_species)} species '
@@ -726,6 +734,9 @@ def determine_species_and_reactions_to_calculate(self) -> bool:
726734
# 1.2. SA
727735
if sa_observables_exist:
728736
species_keys.extend(self.determine_species_based_on_sa())
737+
if self.t3['sensitivity']['global_observables'].lower() == 'idt':
738+
species_idt_keys, rxn_idt_keys = self.determine_params_based_on_sa_idt()
739+
species_keys.extend(species_idt_keys)
729740
# 1.3. collision violators
730741
if self.t3['options']['collision_violators_thermo']:
731742
species_keys.extend(coll_vio_spc_keys)
@@ -734,6 +745,10 @@ def determine_species_and_reactions_to_calculate(self) -> bool:
734745
# 2.1. SA
735746
if sa_observables_exist:
736747
reaction_keys.extend(self.determine_reactions_based_on_sa())
748+
if self.t3['sensitivity']['global_observables'].lower() == 'idt':
749+
if rxn_idt_keys is None:
750+
species_idt_keys, rxn_idt_keys = self.determine_params_based_on_sa_idt()
751+
reaction_keys.extend(rxn_idt_keys)
737752
# 2.2. collision violators
738753
if self.t3['options']['collision_violators_rates']:
739754
reaction_keys.extend(coll_vio_rxn_keys)
@@ -812,6 +827,58 @@ def determine_species_based_on_sa(self) -> List[int]:
812827

813828
return species_keys
814829

830+
def determine_params_based_on_sa_idt(self) -> Tuple[List[int], List[int]]:
831+
"""
832+
Determine species or reactions to calculate based on sensitivity analysis of IDT.
833+
834+
Returns:
835+
List[int]: Entries are T3 species indices of species determined to be calculated based on IDT SA.
836+
"""
837+
visited_species, species_keys = list(), list()
838+
visited_rxns, rxn_keys = list(), list()
839+
if self.sa_dict_idt is None:
840+
self.logger.error(f"T3's sa_dict_idt was None. Please check that the input file contains a proper "
841+
f"'sensitivity' block, that 'IDT' was defined in the global_observables list, "
842+
f"and/or that SA was run successfully.\n"
843+
f"Not performing refinement based on IDT sensitivity analysis!")
844+
return species_keys, rxn_keys
845+
for token in ['thermo', 'kinetics']:
846+
for r, reactor_idt_data in self.sa_dict_idt[token]['IDT'].items():
847+
for phi, phi_data in reactor_idt_data.items():
848+
for p, p_data in phi_data.items():
849+
for t, idt_data in p_data.items():
850+
for index, sa in idt_data.items():
851+
if token == 'thermo':
852+
if index not in visited_species:
853+
visited_species.append(index)
854+
species = self.get_species_by_index(index)
855+
if self.species_requires_refinement(species=species):
856+
reason = f'(i {self.iteration}) IDT is sensitive to this species.'
857+
key = self.add_species(species=species, reasons=reason)
858+
if key is not None:
859+
species_keys.append(key)
860+
elif token == 'kinetics':
861+
if index not in visited_rxns:
862+
visited_rxns.append(index)
863+
reaction = self.get_reaction_by_index(index)
864+
if self.reaction_requires_refinement(reaction):
865+
reason = f'(i {self.iteration}) IDT is sensitive to this reaction.'
866+
key = self.add_reaction(reaction=reaction, reasons=reason)
867+
if key is not None:
868+
rxn_keys.append(key)
869+
for spc in reaction.reactants + reaction.products:
870+
spc_index = self.get_species_key(spc)
871+
if spc_index not in visited_species:
872+
visited_species.append(spc_index)
873+
if self.species_requires_refinement(species=spc):
874+
reason = f'(i {self.iteration}) IDT is sensitive to a reaction ' \
875+
f'({reaction}) in which this species participates.'
876+
key = self.add_species(species=spc, reasons=reason)
877+
if key is not None:
878+
species_keys.append(key)
879+
return species_keys, rxn_keys
880+
881+
815882
def determine_reactions_based_on_sa(self) -> List[int]:
816883
"""
817884
Determine reaction rate coefficients to calculate based on sensitivity analysis.
@@ -1178,6 +1245,21 @@ def get_species_key(self,
11781245
return key
11791246
return None
11801247

1248+
def get_species_by_index(self, index) -> Optional[Species]:
1249+
"""
1250+
Get a species object by its T3 index.
1251+
1252+
Args:
1253+
index (int): The species index.
1254+
1255+
Returns:
1256+
Optional[Species]: The species object if it exists, ``None`` if it does not.
1257+
"""
1258+
for key, species_dict in self.species.items():
1259+
if key == index:
1260+
return species_dict['object']
1261+
return None
1262+
11811263
def get_reaction_key(self,
11821264
reaction: Optional[Reaction] = None,
11831265
label: Optional[str] = None,
@@ -1206,6 +1288,21 @@ def get_reaction_key(self,
12061288
return key
12071289
return None
12081290

1291+
def get_reaction_by_index(self, index) -> Optional[Reaction]:
1292+
"""
1293+
Get a reaction object by its T3 index.
1294+
1295+
Args:
1296+
index (int): The reaction index.
1297+
1298+
Returns:
1299+
Optional[Reaction]: The reaction object if it exists, ``None`` if it does not.
1300+
"""
1301+
for key, reaction_dict in self.reactions.items():
1302+
if key == index:
1303+
return reaction_dict['object']
1304+
return None
1305+
12091306
def load_species_and_reactions_from_chemkin_file(self) -> Tuple[List[Species], List[Reaction]]:
12101307
"""
12111308
Load RMG Species and Reaction objects from the annotated Chemkin file.
@@ -1251,7 +1348,7 @@ def add_species(self,
12511348
) -> Optional[int]:
12521349
"""
12531350
Add a species to self.species and to self.qm['species'].
1254-
If the species already exists in self.species, only the reasons to compute will be appended.
1351+
If the species already exists in self.species, only the reasons to compute it will be appended.
12551352
12561353
Args:
12571354
species (Species): The species to consider.

‎t3/simulate/adapter.py‎

Lines changed: 7 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -5,6 +5,7 @@
55
"""
66

77
from abc import ABC, abstractmethod
8+
from typing import Optional
89

910

1011
class SimulateAdapter(ABC):
@@ -27,7 +28,12 @@ def simulate(self):
2728
pass
2829

2930
@abstractmethod
30-
def get_sa_coefficients(self):
31+
def get_sa_coefficients(self,
32+
top_SA_species: int = 10,
33+
top_SA_reactions: int = 10,
34+
max_workers: int = 24,
35+
save_yaml: bool = True,
36+
) -> Optional[dict]:
3137
"""
3238
Obtain the sensitivity analysis coefficients.
3339

‎t3/simulate/cantera_constantHP.py‎

Lines changed: 14 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -382,12 +382,23 @@ def simulate(self):
382382

383383
self.all_data.append((time, condition_data, reaction_sensitivity_data, thermodynamic_sensitivity_data))
384384

385-
def get_sa_coefficients(self):
385+
def get_sa_coefficients(self,
386+
top_SA_species: int = 10,
387+
top_SA_reactions: int = 10,
388+
max_workers: int = 24,
389+
save_yaml: bool = True,
390+
) -> Optional[dict]:
386391
"""
387392
Obtain the SA coefficients.
388393
394+
Args:
395+
top_SA_species (int, optional): The number of top sensitive species to return.
396+
top_SA_reactions (int, optional): The number of top sensitive reactions to return.
397+
max_workers (int, optional): The maximal number of workers to use for parallel processing.
398+
save_yaml (bool, optional): Save the SA dictionary to a YAML file.
399+
389400
Returns:
390-
sa_dict (dict): a SA dictionary, whose structure is given in the docstring for T3/t3/main.py
401+
sa_dict (Optional[dict]): a SA dictionary, whose structure is given in the docstring for T3/t3/main.py
391402
"""
392403
sa_dict = {'kinetics': dict(), 'thermo': dict(), 'time': list()}
393404

@@ -411,7 +422,7 @@ def get_sa_coefficients(self):
411422
sa_dict['thermo'][observable_label] = dict()
412423
parameter = get_parameter_from_header(spc)
413424
sa_dict['thermo'][observable_label][parameter] = spc.data
414-
425+
# todo: add save yaml to all adapters, use top SA species and reactions, change usage in main
415426
return sa_dict
416427

417428
def get_idt_by_T(self):

‎t3/simulate/cantera_constantTP.py‎

Lines changed: 13 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -381,12 +381,23 @@ def simulate(self):
381381

382382
self.all_data.append((time, condition_data, reaction_sensitivity_data, thermodynamic_sensitivity_data))
383383

384-
def get_sa_coefficients(self):
384+
def get_sa_coefficients(self,
385+
top_SA_species: int = 10,
386+
top_SA_reactions: int = 10,
387+
max_workers: int = 24,
388+
save_yaml: bool = True,
389+
) -> Optional[dict]:
385390
"""
386391
Obtain the SA coefficients.
387392
393+
Args:
394+
top_SA_species (int, optional): The number of top sensitive species to return.
395+
top_SA_reactions (int, optional): The number of top sensitive reactions to return.
396+
max_workers (int, optional): The maximal number of workers to use for parallel processing.
397+
save_yaml (bool, optional): Save the SA dictionary to a YAML file.
398+
388399
Returns:
389-
sa_dict (dict): a SA dictionary, whose structure is given in the docstring for T3/t3/main.py
400+
sa_dict (Optional[dict]): a SA dictionary, whose structure is given in the docstring for T3/t3/main.py
390401
"""
391402
sa_dict = {'kinetics': dict(), 'thermo': dict(), 'time': list()}
392403

‎t3/simulate/cantera_constantUV.py‎

Lines changed: 13 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -381,12 +381,23 @@ def simulate(self):
381381

382382
self.all_data.append((time, condition_data, reaction_sensitivity_data, thermodynamic_sensitivity_data))
383383

384-
def get_sa_coefficients(self):
384+
def get_sa_coefficients(self,
385+
top_SA_species: int = 10,
386+
top_SA_reactions: int = 10,
387+
max_workers: int = 24,
388+
save_yaml: bool = True,
389+
) -> Optional[dict]:
385390
"""
386391
Obtain the SA coefficients.
387392
393+
Args:
394+
top_SA_species (int, optional): The number of top sensitive species to return.
395+
top_SA_reactions (int, optional): The number of top sensitive reactions to return.
396+
max_workers (int, optional): The maximal number of workers to use for parallel processing.
397+
save_yaml (bool, optional): Save the SA dictionary to a YAML file.
398+
388399
Returns:
389-
sa_dict (dict): a SA dictionary, whose structure is given in the docstring for T3/t3/main.py
400+
sa_dict (Optional[dict]): a SA dictionary, whose structure is given in the docstring for T3/t3/main.py
390401
"""
391402
sa_dict = {'kinetics': dict(), 'thermo': dict(), 'time': list()}
392403

‎t3/simulate/rmg_constantTP.py‎

Lines changed: 16 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -18,6 +18,8 @@
1818
from rmgpy.tools.loader import load_rmg_py_job
1919
from rmgpy.tools.plot import plot_sensitivity
2020

21+
from arc.common import read_yaml_file, save_yaml_file
22+
2123
from t3.common import get_chem_to_rmg_rxn_index_map, get_species_by_label, get_values_within_range, \
2224
get_observable_label_from_header, get_parameter_from_header, time_lapse
2325
from t3.simulate.adapter import SimulateAdapter
@@ -224,10 +226,21 @@ def simulate(self):
224226

225227
self.logger.info(f'Simulation via RMG completed, execution time: {time_lapse(tic)}')
226228

227-
def get_sa_coefficients(self) -> Optional[dict]:
229+
def get_sa_coefficients(self,
230+
top_SA_species: int = 10,
231+
top_SA_reactions: int = 10,
232+
max_workers: int = 24,
233+
save_yaml: bool = True,
234+
) -> Optional[dict]:
228235
"""
229236
Obtain the SA coefficients.
230237
238+
Args:
239+
top_SA_species (int, optional): The number of top sensitive species to return.
240+
top_SA_reactions (int, optional): The number of top sensitive reactions to return.
241+
max_workers (int, optional): The maximal number of workers to use for parallel processing.
242+
save_yaml (bool, optional): Save the SA dictionary to a YAML file.
243+
231244
Returns:
232245
sa_dict (Optional[dict]): An SA dictionary, structure is given in the docstring for T3/t3/main.py
233246
"""
@@ -265,6 +278,8 @@ def get_sa_coefficients(self) -> Optional[dict]:
265278
parameter = chem_to_rmg_rxn_index_map[int(parameter)] \
266279
if all(c.isdigit() for c in parameter) else parameter
267280
sa_dict[sa_type][observable_label][parameter] = df[header].values
281+
if save_yaml:
282+
save_yaml_file(path=self.paths['SA dict'], content=sa_dict)
268283
return sa_dict
269284

270285
def generate_rmg_reactors_for_simulation(self) -> dict:

0 commit comments

Comments
 (0)