diff --git a/.gitignore b/.gitignore index e76ad0a..d00f7a0 100644 --- a/.gitignore +++ b/.gitignore @@ -70,3 +70,6 @@ target/ # project-specific data files *.h5 + +# IDE specific files +.vscode diff --git a/pyteck/eval_model.py b/pyteck/eval_model.py index b158280..cde8553 100644 --- a/pyteck/eval_model.py +++ b/pyteck/eval_model.py @@ -4,31 +4,219 @@ # Standard libraries import os -from os.path import splitext, basename import multiprocessing import warnings -import numpy +import numpy as np from scipy.interpolate import UnivariateSpline +import yaml -try: - import yaml -except ImportError: - print('Warning: YAML must be installed to read input file.') - raise - -from pyked.chemked import ChemKED, DataPoint +from pyked.chemked import ChemKED, IgnitionDataPoint, SpeciesProfileDataPoint # Local imports from .utils import units -from .simulation import Simulation +from .simulation import AutoIgnitionSimulation, JSRSimulation + -min_deviation = 0.10 +min_deviation = 0.1 """float: minimum allowable standard deviation for experimental data""" -def create_simulations(dataset, properties): - """Set up individual simulations for each ignition delay value. +def ignition_dataset_processing(results, print_results=False): + """Function to process all of the simulated datapoints + from a single ignition delay dataset + + results is a list of pyteck.simulation.AutoIgnitionSimulation + """ + dataset_meta = {} + + ignition_delays_exp = np.zeros(len(results)) + ignition_delays_sim = np.zeros(len(results)) + + ############################################# + # Determine standard deviation of the dataset + ############################################# + ign_delay = [ + sim.properties.ignition_delay.to('second').value.magnitude + if hasattr(sim.properties.ignition_delay, 'value') + else sim.properties.ignition_delay.to('second').magnitude + for sim in results + ] + + datapoints = [sim.properties for sim in results] + + # get variable that is changing across datapoints + # variable = get_changing_variable(properties.datapoints) + variable = get_changing_variable(datapoints) + + # for ignition delay, use logarithm of values + standard_dev = estimate_std_dev(variable, np.log(ign_delay)) + dataset_meta['standard deviation'] = float(standard_dev) + dataset_meta['datapoints'] = [] + # dataset_meta[''] add the dataset name + + for idx, sim in enumerate(results): + + sim.process_results() + dataset_meta.update(sim.meta) + if hasattr(sim.properties.ignition_delay, 'value'): + ignition_delay = sim.properties.ignition_delay.value + else: + ignition_delay = sim.properties.ignition_delay + + if hasattr(ignition_delay, 'nominal_value'): + ignition_delay = ignition_delay.nominal_value * units.second + + dataset_meta['datapoints'].append({ + 'experimental ignition delay': str(ignition_delay), + 'simulated ignition delay': str(sim.meta['simulated-ignition-delay']), + 'temperature': str(sim.properties.temperature), + 'pressure': str(sim.properties.pressure), + 'composition': [{ + 'InChI': sim.properties.composition[spec].InChI, + 'species-name': sim.properties.composition[spec].species_name, + 'amount': str(sim.properties.composition[spec].amount.magnitude), + } for spec in sim.properties.composition], + 'composition type': sim.properties.composition_type, + }) + + ignition_delays_exp[idx] = ignition_delay.magnitude + ignition_delays_sim[idx] = sim.meta['simulated-ignition-delay'].magnitude + + # calculate error function for this dataset + error_func = np.power( + (np.log(ignition_delays_sim) - np.log(ignition_delays_exp)) + / standard_dev, 2 + ) + error_func = np.nanmean(error_func) + dataset_meta['error function'] = float(error_func) + + dev_func = ( + np.log(ignition_delays_sim) + - np.log(ignition_delays_exp) + ) / standard_dev + dev_func = np.nanmean(dev_func) + dataset_meta['absolute deviation'] = float(dev_func) + + return dataset_meta + + +def JSR_dataset_processing(results, print_results=False): + """Function to process all of the simulated datapoints + from a single jet-stirred reactor dataset + + results is a list of pyteck.simulation.JSRSimulation + """ + dataset_meta = {} + dataset_meta['datapoints'] = [] + expt_target_species_profiles = [] + simulated_species_profiles = [] + inlet_temperatures = [] + for i, sim in enumerate(results): + dataset_meta.update(sim.meta) + concentration = sim.process_results() + species_name = sim.meta['species_name'] + + inlet_composition = {} + for k, v in sim.properties.inlet_composition.items(): + inlet_composition[k] = v.amount.magnitude.nominal_value + expt_target_species_profile = sim.properties.outlet_composition[species_name].amount + inlet_temperature = sim.properties.temperature + dataset_meta['datapoints'].append({ + # 'inlet composition': str(inlet_composition), + 'experimental species profile': str(expt_target_species_profile), + 'simulated species profile': str(concentration), + 'temperature': str(sim.properties.temperature), + 'pressure': str(sim.properties.pressure), + }) + + expt_target_species_profiles.append(expt_target_species_profile.magnitude) + simulated_species_profiles.append(concentration) + inlet_temperatures.append(inlet_temperature) + + # calculate error function for this dataset + experimental_trapz = np.trapz(inlet_temperatures, expt_target_species_profiles) + simulated_trapz = np.trapz(inlet_temperatures, simulated_species_profiles) + if print_results: + print("Difference between AUC:{}".format(experimental_trapz - simulated_trapz)) + return dataset_meta + + +def ignition_total_processing(results_stats, print_results=False): + """Function to get overall statistics for ignition delays + from all provided datasets and simulations + """ + output = {'datasets': []} + # NOTE results_stats already excludes skipped datasets + error_func_sets = np.zeros(len(results_stats)) + dev_func_sets = np.zeros(len(results_stats)) + for i, dataset_meta in enumerate(results_stats): + dev_func_sets[i] = dataset_meta['absolute deviation'] + error_func_sets[i] = dataset_meta['error function'] + output['datasets'].append(dataset_meta) + + # Overall error function + error_func = np.nanmean(error_func_sets) + if print_results: + print('overall error function: ' + repr(error_func)) + print('error standard deviation: ' + repr(np.nanstd(error_func_sets))) + + # Absolute deviation function + abs_dev_func = np.nanmean(dev_func_sets) + if print_results: + print('absolute deviation function: ' + repr(abs_dev_func)) + + output['average error function'] = float(error_func) + output['error function standard deviation'] = float(np.nanstd(error_func_sets)) + output['average deviation function'] = float(abs_dev_func) + + return output + + +def JSR_total_processing(results_stats, print_results=False): + """Function to get overall statistics for species concentrations + from all provided datasets and simulations + """ + return results_stats + + +def SimulationFactory(datapoint_class): + """ + Get the Simulation that corresponds to the input datapoint type + """ + simulations = { + IgnitionDataPoint: AutoIgnitionSimulation, + SpeciesProfileDataPoint: JSRSimulation, + } + return simulations[datapoint_class] + + +def DatasetProcessingFactory(datapoint_class): + """ + Get the dataset statistics post-processing function + that corresponds to the input datapoint type + """ + simulations = { + IgnitionDataPoint: ignition_dataset_processing, + SpeciesProfileDataPoint: JSR_dataset_processing, + } + return simulations[datapoint_class] + + +def TotalProcessingFactory(datapoint_class): + """ + Get the total statistics post-processing function that + corresponds to the input datapoint type + """ + simulations = { + IgnitionDataPoint: ignition_total_processing, + SpeciesProfileDataPoint: JSR_total_processing, + } + return simulations[datapoint_class] + + +def create_simulations(dataset, properties, **kwargs): + """Set up individual simulations for each experimental datapoint. Parameters ---------- @@ -40,25 +228,32 @@ def create_simulations(dataset, properties): Returns ------- simulations : list - List of :class:`Simulation` objects for each simulation + List of :class:`AutoignitionSimulation` objects for each simulation or + List of :class:`JSRSimulation` objects for each simulation """ - simulations = [] for idx, case in enumerate(properties.datapoints): sim_meta = {} + sim_meta.update(kwargs) # Common metadata sim_meta['data-file'] = dataset - sim_meta['id'] = splitext(basename(dataset))[0] + '_' + str(idx) - - simulations.append(Simulation(properties.experiment_type, - properties.apparatus.kind, - sim_meta, - case - ) - ) + sim_meta['id'] = os.path.splitext(os.path.basename(dataset))[0] + '_' + str(idx) + Simulation = SimulationFactory(type(case)) + + simulations.append( + Simulation( + properties.experiment_type, + properties.apparatus.kind, + sim_meta, + case, + **kwargs, + ) + ) + return simulations + def simulation_worker(sim_tuple): """Worker for multiprocessing of simulation cases. @@ -70,16 +265,18 @@ def simulation_worker(sim_tuple): Returns ------- - sim : ``Simulation`` - Simulation case with calculated ignition delay. + sim : ``AutoignitionSimulation`` + or ``JSRSimulation`` + Simulation case with calculated results """ sim, model_file, model_spec_key, path, restart = sim_tuple + simulation_type = type(sim) sim.setup_case(model_file, model_spec_key, path) sim.run_case(restart) - sim = Simulation(sim.kind, sim.apparatus, sim.meta, sim.properties) + sim = simulation_type(sim.kind, sim.apparatus, sim.meta, sim.properties) return sim @@ -99,20 +296,18 @@ def estimate_std_dev(indep_variable, dep_variable): Standard deviation of difference between data and best-fit line """ - assert len(indep_variable) == len(dep_variable), \ 'independent and dependent variables not the same length' # ensure no repetition of independent variable by taking average of associated dependent # variables and removing duplicates - vals, count = numpy.unique(indep_variable, return_counts=True) + vals, count = np.unique(indep_variable, return_counts=True) repeated = vals[count > 1] for val in repeated: - idx, = numpy.where(indep_variable == val) - dep_variable[idx[0]] = numpy.mean(dep_variable[idx]) - dep_variable = numpy.delete(dep_variable, idx[1:]) - indep_variable = numpy.delete(indep_variable, idx[1:]) - + idx, = np.where(indep_variable == val) + dep_variable[idx[0]] = np.mean(dep_variable[idx]) + dep_variable = np.delete(dep_variable, idx[1:]) + indep_variable = np.delete(indep_variable, idx[1:]) # ensure data sorted based on independent variable to avoid some problems sorted_vars = sorted(zip(indep_variable, dep_variable)) @@ -128,7 +323,7 @@ def estimate_std_dev(indep_variable, dep_variable): else: spline = UnivariateSpline(indep_variable, dep_variable) - standard_dev = numpy.std(dep_variable - spline(indep_variable)) + standard_dev = np.std(dep_variable - spline(indep_variable)) if standard_dev < min_deviation: print('Standard deviation of {:.2f} too low, ' @@ -143,8 +338,8 @@ def get_changing_variable(cases): Parameters ---------- - cases : list(pyked.chemked.DataPoint) - List of DataPoint with experimental case data. + cases : list(pyked.chemked.IgnitionDataPoint) + List of IgnitionDataPoint with experimental case data. Returns ------- @@ -176,24 +371,36 @@ def get_changing_variable(cases): changing_var = 'temperature' if changing_var == 'temperature': - variable = [case.temperature.value.magnitude if hasattr(case.temperature, 'value') - else case.temperature.magnitude - for case in cases - ] + variable = [ + case.temperature.value.magnitude if hasattr(case.temperature, 'value') + else case.temperature.magnitude + for case in cases + ] elif changing_var == 'pressure': - variable = [case.pressure.value.magnitude if hasattr(case.pressure, 'value') - else case.pressure.magnitude - for case in cases - ] + variable = [ + case.pressure.value.magnitude if hasattr(case.pressure, 'value') + else case.pressure.magnitude + for case in cases + ] + if variable[0].__class__.__name__ == 'Variable': + variable = [var.nominal_value for var in variable] return variable -def evaluate_model(model_name, spec_keys_file, dataset_file, - data_path='data', model_path='models', - results_path='results', model_variant_file=None, - num_threads=None, print_results=False, restart=False, - skip_validation=False, - ): +def evaluate_model( + model_name, + spec_keys_file, + dataset_file, + data_path='data', + model_path='models', + results_path='results', + model_variant_file=None, + species_name=None, + num_threads=None, + print_results=False, + restart=False, + skip_validation=False +): """Evaluates the ignition delay error of a model for a given dataset. Parameters @@ -244,72 +451,51 @@ def evaluate_model(model_name, spec_keys_file, dataset_file, with open(model_variant_file, 'r') as f: model_variant = yaml.safe_load(f) - # Read dataset list + # Read dataset list, ignoring blank lines and lines starting with # + dataset_list = [] with open(dataset_file, 'r') as f: - dataset_list = f.read().splitlines() - - error_func_sets = numpy.zeros(len(dataset_list)) - dev_func_sets = numpy.zeros(len(dataset_list)) - - # Dictionary with all output data - output = {'model': model_name, 'datasets': []} + lines = f.readlines() + for line in lines: + formatted_line = line.strip() + if formatted_line == '' or formatted_line[0] == '#': + continue + dataset_list.append(formatted_line) # If number of threads not specified, use either max number of available # cores minus 1, or use 1 if multiple cores not available. if not num_threads: - num_threads = multiprocessing.cpu_count()-1 or 1 + num_threads = multiprocessing.cpu_count() - 1 or 1 # Loop through all datasets + skipped_datasets = [] + results_stats = [] for idx_set, dataset in enumerate(dataset_list): - dataset_meta = {'dataset': dataset, 'dataset_id': idx_set} + dataset_meta = { + 'model_name': model_name, + 'species_name': species_name, + } # Create individual simulation cases for each datapoint in this set properties = ChemKED(os.path.join(data_path, dataset), skip_validation=skip_validation) - simulations = create_simulations(dataset, properties) - - ignition_delays_exp = numpy.zeros(len(simulations)) - ignition_delays_sim = numpy.zeros(len(simulations)) - - ############################################# - # Determine standard deviation of the dataset - ############################################# - ign_delay = [case.ignition_delay.to('second').value.magnitude - if hasattr(case.ignition_delay, 'value') - else case.ignition_delay.to('second').magnitude - for case in properties.datapoints - ] - - # get variable that is changing across datapoints - variable = get_changing_variable(properties.datapoints) - # for ignition delay, use logarithm of values - standard_dev = estimate_std_dev(variable, numpy.log(ign_delay)) - dataset_meta['standard deviation'] = float(standard_dev) ####################################################### # Need to check if Ar or He in reactants but not model, # and if so skip this dataset (for now). ####################################################### - if ((any(['Ar' in spec for case in properties.datapoints - for spec in case.composition] - ) - and 'Ar' not in model_spec_key[model_name] - ) or - (any(['He' in spec for case in properties.datapoints - for spec in case.composition] - ) - and 'He' not in model_spec_key[model_name] - ) - ): - warnings.warn('Warning: Ar or He in dataset, but not in model. Skipping.', - RuntimeWarning - ) - error_func_sets[idx_set] = numpy.nan + Ar_in_model = 'Ar' in model_spec_key[model_name] + He_in_model = 'He' in model_spec_key[model_name] + Ar_in_dataset = any(case.species_in_datapoint('Ar') for case in properties.datapoints) + He_in_dataset = any(case.species_in_datapoint('He') for case in properties.datapoints) + if (Ar_in_dataset and not Ar_in_model) or (He_in_dataset and not He_in_model): + warnings.warn( + 'Warning: Ar or He in dataset, but not in model. Skipping.', + RuntimeWarning + ) + skipped_datasets.append(idx_set) # TODO set the error to NAN continue - # Use available number of processors minus one, - # or one process if single core. - pool = multiprocessing.Pool(processes=num_threads) + simulations = create_simulations(dataset, properties, **dataset_meta) # setup all cases jobs = [] @@ -323,7 +509,7 @@ def evaluate_model(model_name, spec_keys_file, dataset_file, bath_gases = set(model_variant[model_name]['bath gases']) gases = bath_gases.intersection( set([c['species-name'] for c in sim.properties.composition]) - ) + ) # If only one bath gas present, use that. If multiple, use the # predominant species. If none of the designated bath gases @@ -347,11 +533,10 @@ def evaluate_model(model_name, spec_keys_file, dataset_file, # choose closest pressure # better way to do this? - i = numpy.argmin(numpy.abs(numpy.array( - [float(n) - for n in list(model_variant[model_name]['pressures']) - ] - ) - pres)) + i = np.argmin(np.abs(np.array([ + float(n) + for n in list(model_variant[model_name]['pressures']) + ]) - pres)) pres = list(model_variant[model_name]['pressures'])[i] model_mod += model_variant[model_name]['pressures'][pres] @@ -361,80 +546,35 @@ def evaluate_model(model_name, spec_keys_file, dataset_file, jobs.append([sim, model_file, model_spec_key[model_name], results_path, restart]) - # run all cases - jobs = tuple(jobs) - results = pool.map(simulation_worker, jobs) - - # not adding more proceses, and ensure all finished - pool.close() - pool.join() - - dataset_meta['datapoints'] = [] - - for idx, sim in enumerate(results): - sim.process_results() - - if hasattr(sim.properties.ignition_delay, 'value'): - ignition_delay = sim.properties.ignition_delay.value - else: - ignition_delay = sim.properties.ignition_delay - - if hasattr(ignition_delay, 'nominal_value'): - ignition_delay = ignition_delay.nominal_value * units.second - - dataset_meta['datapoints'].append( - {'experimental ignition delay': str(ignition_delay), - 'simulated ignition delay': str(sim.meta['simulated-ignition-delay']), - 'temperature': str(sim.properties.temperature), - 'pressure': str(sim.properties.pressure), - 'composition': [{'InChI': sim.properties.composition[spec].InChI, - 'species-name': sim.properties.composition[spec].species_name, - 'amount': str(sim.properties.composition[spec].amount.magnitude), - } for spec in sim.properties.composition], - 'composition type': sim.properties.composition_type, - }) - - ignition_delays_exp[idx] = ignition_delay.magnitude - ignition_delays_sim[idx] = sim.meta['simulated-ignition-delay'].magnitude - - # calculate error function for this dataset - error_func = numpy.power( - (numpy.log(ignition_delays_sim) - - numpy.log(ignition_delays_exp)) / standard_dev, 2 - ) - error_func = numpy.nanmean(error_func) - error_func_sets[idx_set] = error_func - dataset_meta['error function'] = float(error_func) - - dev_func = (numpy.log(ignition_delays_sim) - - numpy.log(ignition_delays_exp) - ) / standard_dev - dev_func = numpy.nanmean(dev_func) - dev_func_sets[idx_set] = dev_func - dataset_meta['absolute deviation'] = float(dev_func) - - output['datasets'].append(dataset_meta) - + if num_threads == 1: + # Don't use the threadpool if only 1 processor (useful for debugging) + results = [] + for job in jobs: + results.append(simulation_worker(job)) + else: + pool = multiprocessing.Pool(processes=num_threads) + jobs = tuple(jobs) + results = pool.map(simulation_worker, jobs) + + # not adding more proceses, and ensure all finished + pool.close() + pool.join() + + # process the results for each dataset + post_proc = DatasetProcessingFactory(type(properties.datapoints[0])) + results_stats.append(post_proc(results)) if print_results: print('Done with ' + dataset) - # Overall error function - error_func = numpy.nanmean(error_func_sets) - if print_results: - print('overall error function: ' + repr(error_func)) - print('error standard deviation: ' + repr(numpy.nanstd(error_func_sets))) - - # Absolute deviation function - abs_dev_func = numpy.nanmean(dev_func_sets) - if print_results: - print('absolute deviation function: ' + repr(abs_dev_func)) - - output['average error function'] = float(error_func) - output['error function standard deviation'] = float(numpy.nanstd(error_func_sets)) - output['average deviation function'] = float(abs_dev_func) + # process the total stats for multiple datasets + total_proc = TotalProcessingFactory(type(properties.datapoints[0])) + output = total_proc(results_stats) # Write data to YAML file - with open(splitext(basename(model_name))[0] + '-results.yaml', 'w') as f: + with open(os.path.join( + results_path, + os.path.splitext(os.path.basename(model_name))[0] + '-results.yaml' + ), 'w') as f: yaml.dump(output, f) return output diff --git a/pyteck/plotting.py b/pyteck/plotting.py new file mode 100644 index 0000000..e68c81d --- /dev/null +++ b/pyteck/plotting.py @@ -0,0 +1,121 @@ +import pandas as pd +import h5py +import matplotlib.pyplot as plt +import yaml +import cantera as ct +import os +import glob + +from matplotlib.backends.backend_pdf import PdfPages + + +def generate_plots_jsr(model_name, model_path, results_path, spec_keys_file, data_path, plot_path): + + """Generates plots of steady-state concentration over temperature for each species in the species key. + + Parameters + ---------- + model_name : str + Chemical kinetic model filename + model_path : str + Local path for model file. Optional; default = 'models' + results_path : str + Local path for creating results files. Optional; default = 'results' + spec_keys_file : str + Name of YAML file identifying important species + data_path : str + Local path for data files. Optional; default = 'data' + plot_path : bool + Local path for creating the plots pdf. Optional; default = 'jsr_plots' + + """ + + sol = ct.Solution(os.path.join(model_path, model_name)) + h5_list = glob.glob(os.path.join(results_path, '*.h5')) + experimental_files = os.listdir(data_path) + + spec_names_model = (sol.species_names) + + for file in experimental_files: + if os.path.splitext(file)[-1] == ".yaml": # Any none .yaml files are skipped over + if os.path.exists(os.path.join(data_path, file)): + print(f"Loading {file}") + with open(os.path.join(data_path, file), 'r') as f: + data = yaml.load(f, Loader=yaml.SafeLoader) + + else: + raise OSError(f"Couldn't find {os.path.join(data_path,file)}") + else: + print(f"Ignoring none .yaml file {file}") + continue + + if os.path.exists(spec_keys_file): + with open(spec_keys_file, 'r') as k: + key = yaml.load(k, Loader=yaml.SafeLoader) + else: + raise OSError(f"Couldn't find {spec_keys_file}") + + species = data['datapoints'][0]['outlet-composition']['species'] + # experimental file should contain the path to the corresponding csv file + csvfile = data['datapoints'][0]['csvfile'] + + if os.path.exists(csvfile): + exp = pd.read_csv(csvfile) + print(f"Loading {csvfile}") + elif os.path.exists(os.path.join('data', csvfile)): + exp = pd.read_csv(os.path.join('data', csvfile)) + print(f"Loading {os.path.join('data', csvfile)}") + elif os.path.exists(os.path.join(data_path, csvfile)): + exp = pd.read_csv(os.path.join(data_path, csvfile)) + print(f"Loading {os.path.join(data_path, csvfile)}") + else: + raise OSError(f"Couldn't find {csvfile}") + + with PdfPages(os.path.join(plot_path, 'jsr_plots') + '-' + model_name + '.pdf') as plot_pdf: + for sp in species: + name_in_data = sp['species-name'] + + if name_in_data in key[model_name].keys(): + + name_in_model = key[model_name][sp['species-name']] + print('Plotting concentration for ' + name_in_model) + + temps = [] + concs = [] + + # get the position in which this species is listed in the model + i = 0 + for name in spec_names_model: + if(name == name_in_model): + break + else: + i = i + 1 + # iterate through all results files (each for a single temp) + for h5 in h5_list: + f = h5py.File(os.path.join(results_path, h5), 'r') + dset = f['simulation'] + + temp = dset[1][1] + conc = (dset[-1][4][i]) + + concs.append(conc) + temps.append(temp) + + f.close() + + temps, concs = zip(*sorted(zip(temps, concs))) # sort both lists + + plt.figure() + plt.plot(temps, concs, linestyle='solid') + plt.title(sp['species-name'] + ' concentration') + + temps = exp['Temperature'] + concs = exp[sp['species-name']] + + plt.scatter(temps, concs) + plt.legend(['simulated', 'experimental']) + plt.xlabel('Temperature (K)') + plt.ylabel('Mole Fraction') + + plot_pdf.savefig() + plt.close() diff --git a/pyteck/simulation.py b/pyteck/simulation.py index ec37087..188a2da 100644 --- a/pyteck/simulation.py +++ b/pyteck/simulation.py @@ -1,6 +1,9 @@ """ -.. moduleauthor:: Kyle Niemeyer +.. moduleauthor:: Kyle Niemeyer , + Sai Krishna Sirumalla + Sevy Harris + """ # Python 2 compatibility @@ -9,10 +12,10 @@ # Standard libraries import os -from collections import namedtuple import warnings import numpy as np -from datetime import datetime + +from abc import ABC, abstractmethod # Related modules try: @@ -32,6 +35,7 @@ from .utils import units from .detect_peaks import detect_peaks + def first_derivative(x, y): """Evaluates first derivative using second-order finite differences. @@ -39,11 +43,11 @@ def first_derivative(x, y): one-sided difference at boundaries. :param x: Independent variable array - :type x: np.ndarray + :type x: numpy.ndarray :param y: Dependent variable array - :type y: np.ndarray + :type y: numpy.ndarray :return: First derivative, :math:`dy/dx` - :rtype: np.ndarray + :rtype: numpy.ndarray """ return np.gradient(y, x, edge_order=2) @@ -56,7 +60,7 @@ def sample_rising_pressure(time_end, init_pres, freq, pressure_rise_rate): :param float freq: Frequency of sampling, in Hz :param float pressure_rise_rate: Pressure rise rate, in s^-1 :return: List of times and pressures - :rtype: list of np.ndarray + :rtype: list of numpy.ndarray """ times = np.arange(0.0, time_end + (1.0 / freq), (1.0 / freq)) pressures = init_pres * (pressure_rise_rate * times + 1.0) @@ -73,7 +77,7 @@ def create_volume_history(mech, temp, pres, reactants, pres_rise, time_end): :param float pres_rise: Pressure rise rate, in s^-1 :param float time_end: End time of simulation in s :return: List of times and volumes - :rtype: list of np.ndarray + :rtype: list of numpy.ndarray """ gas = ct.Solution(mech) gas.TPX = temp, pres, reactants @@ -93,77 +97,6 @@ def create_volume_history(mech, temp, pres, reactants, pres_rise, time_end): return [times, volumes] -def get_ignition_delay(time, target, target_name, ignition_type): - """Identify ignition delay based on time, target, and type of detection. - """ - if ignition_type == 'max': - # Get indices of peaks - peak_inds = detect_peaks(target, edge=None, mph=1.e-9*np.max(target)) - - # Get index of largest peak (overall ignition delay) - max_ind = peak_inds[np.argmax(target[peak_inds])] - - #ign_delays = time[peak_inds[np.where((time[peak_inds[peak_inds <= max_ind]]) > 0.0)]] - - ign_delays = time[peak_inds[peak_inds <= max_ind]] - - elif ignition_type == 'd/dt max': - target = first_derivative(time, target) - # Get indices of peaks. Set a minimum peak height of 1e-7% of the - # maximum value to avoid noise peaks. - peak_inds = detect_peaks(target, edge=None, mph=1.e-9*np.max(target)) - - # Get index of largest peak (overall ignition delay) - max_ind = peak_inds[np.argmax(target[peak_inds])] - - ign_delays = time[peak_inds[np.where((time[peak_inds[peak_inds <= max_ind]]) > 0.0)]] - - elif ignition_type == '1/2 max': - # maximum value, and associated index - max_val = np.max(target) - peak_inds = detect_peaks(target, edge=None, mph=1.e-9*np.max(target)) - max_ind = peak_inds[np.argmax(target[peak_inds])] - - # TODO: interpolate for actual half-max value - # Find index associated with the 1/2 max value, but only consider - # points before the peak - half_idx = (np.abs(target[0:max_ind] - 0.5 * max_val)).argmin() - ign_delays = np.array([time[half_idx]]) - - elif ignition_type == 'd/dt max extrapolated': - # First need to evaluate derivative of the target - target_derivative = first_derivative(time, target) - - # Get indices of peaks, and index of largest peak, which corresponds to - # the point of maximum deriative - peak_inds = detect_peaks(target_derivative, edge=None, mph=1.e-9*np.max(target)) - max_ind = peak_inds[np.argmax(target_derivative[peak_inds])] - - # use slope to extrapolate to intercept with baseline value (0 by default) - ign_delays = np.array([time[max_ind] - (target[max_ind] / target_derivative[max_ind])]) - - # TODO: handle target with nonzero baseline? - else: - warnings.warn('Unable to process ignition type ' + ignition_type + - ', setting result to 0 and continuing', RuntimeWarning - ) - return np.array([0.0]) - - - # something has gone wrong if there is still no peak. This shouldn't be necessary. - if ign_delays.size == 0: - filename = 'target-data-' + str(datetime.now().strftime('%Y_%m_%d_%H_%M_%S')) + '.out' - warnings.warn('No peak found, dumping target data to ' + - filename + ' and continuing', RuntimeWarning - ) - np.savetxt(filename, np.c_[time, target], - header=('time, target (' + target_name + ')') - ) - return np.array([0.0]) - - return ign_delays - - class VolumeProfile(object): """Set the velocity of reactor moving wall via specified volume profile. @@ -188,13 +121,14 @@ def __init__(self, volume_history): :param VolumeHistory volume_history: time and volume history """ - # The time and volume are each stored as a ``np.array`` in the + # The time and volume are each stored as a ``numpy.array`` in the # properties dictionary. The volume is normalized by the first volume # element so that a unit area can be used to calculate the velocity. self.times = volume_history.time.magnitude - volumes = (volume_history.quantity.magnitude / - volume_history.quantity.magnitude[0] - ) + volumes = ( + volume_history.quantity.magnitude + / volume_history.quantity.magnitude[0] + ) # The velocity is calculated by the second-order central differences. self.velocity = first_derivative(self.times, volumes) @@ -251,23 +185,25 @@ def __init__(self, mech_filename, initial_temp, initial_pres, """ [self.times, volumes] = create_volume_history( - mech_filename, initial_temp, initial_pres, - reactants, pressure_rise, time_end - ) + mech_filename, initial_temp, initial_pres, + reactants, pressure_rise, time_end + ) # Calculate velocity by second-order finite difference self.velocity = first_derivative(self.times, volumes) -class Simulation(object): - """Class for ignition delay simulations.""" - - def __init__(self, kind, apparatus, meta, properties): +class Simulation(ABC): + """ + Superclass for simulations. Each simulation must implement + its own setup_case, run_case, and process_results methods + """ + def __init__(self, kind, apparatus, meta, properties, **kwargs): """Initialize simulation case. - :param kind: Kind of experiment (e.g., 'ignition delay') + :param kind: Kind of experiment (e.g., 'ignition delay' or 'species profile') :type kind: str - :param apparatus: Type of apparatus ('shock tube' or 'rapid compression machine') + :param apparatus: Type of apparatus ('shock tube','rapid compression machine' or jet-stirred reactor ') :type apparatus: str :param meta: some metadata for this case :type meta: dict @@ -279,6 +215,26 @@ def __init__(self, kind, apparatus, meta, properties): self.meta = meta self.properties = properties + @abstractmethod + def setup_case(self): + raise NotImplementedError + + @abstractmethod + def run_case(self, restart=False): + raise NotImplementedError + + @abstractmethod + def process_results(self): + raise NotImplementedError + + +class AutoIgnitionSimulation(Simulation): + """Class for shock tube and rapid compression machine + simulations with the goal of computing ignition delay times + """ + def __init__(self, *args, **kwargs): + super(self.__class__, self).__init__(*args, **kwargs) + def setup_case(self, model_file, species_key, path=''): """Sets up the simulation case to be run. @@ -305,10 +261,11 @@ def setup_case(self, model_file, species_key, path=''): self.properties.pressure.ito('pascal') # convert reactant names to those needed for model - reactants = [species_key[self.properties.composition[spec].species_name] + ':' + - str(self.properties.composition[spec].amount.magnitude) - for spec in self.properties.composition - ] + reactants = [ + species_key[self.properties.composition[spec].species_name] + ':' + + str(self.properties.composition[spec].amount.magnitude) + for spec in self.properties.composition + ] reactants = ','.join(reactants) # need to extract values from Quantity or Measurement object @@ -318,6 +275,7 @@ def setup_case(self, model_file, species_key, path=''): temp = self.properties.temperature.nominal_value else: temp = self.properties.temperature.magnitude + if hasattr(self.properties.pressure, 'value'): pres = self.properties.pressure.value.magnitude elif hasattr(self.properties.pressure, 'nominal_value'): @@ -332,7 +290,6 @@ def setup_case(self, model_file, species_key, path=''): self.gas.TPY = temp, pres, reactants else: raise(BaseException('error: not supported')) - return # Create non-interacting ``Reservoir`` on other side of ``Wall`` env = ct.Reservoir(ct.Solution('air.xml')) @@ -353,34 +310,38 @@ def setup_case(self, model_file, species_key, path=''): else: pres_rise = self.properties.pressure_rise.magnitude - self.wall = ct.Wall(self.reac, env, A=1.0, - velocity=PressureRiseProfile( - model_file, - self.gas.T, - self.gas.P, - self.gas.X, - pres_rise, - self.time_end - ) - ) - - elif (self.apparatus == 'rapid compression machine' and - self.properties.volume_history is None - ): + self.wall = ct.Wall( + self.reac, env, A=1.0, + velocity=PressureRiseProfile( + model_file, + self.gas.T, + self.gas.P, + self.gas.X, + pres_rise, + self.time_end + ) + ) + + elif ( + self.apparatus == 'rapid compression machine' + and self.properties.volume_history is None + ): # Rapid compression machine modeled by constant UV self.wall = ct.Wall(self.reac, env, A=1.0, velocity=0) - elif (self.apparatus == 'rapid compression machine' and - self.properties.volume_history is not None - ): + elif ( + self.apparatus == 'rapid compression machine' + and self.properties.volume_history is not None + ): # Rapid compression machine modeled with volume-time history # First convert time units if necessary self.properties.volume_history.time.ito('second') - self.wall = ct.Wall(self.reac, env, A=1.0, - velocity=VolumeProfile(self.properties.volume_history) - ) + self.wall = ct.Wall( + self.reac, env, A=1.0, + velocity=VolumeProfile(self.properties.volume_history) + ) # Number of solution variables is number of species + mass, # volume, temperature @@ -393,7 +354,7 @@ def setup_case(self, model_file, species_key, path=''): if self.properties.volume_history is not None: # Minimum difference between volume profile times min_time = np.min(np.diff(self.properties.volume_history.time.magnitude)) - self.reac_net.set_max_time_step(min_time) + self.reac_net.max_time_step = min_time # Check if species ignition target, that species is present. if self.properties.ignition_type['target'] not in ['pressure', 'temperature']: @@ -423,7 +384,7 @@ def setup_case(self, model_file, species_key, path=''): warnings.warn( spec + ' not found in model; falling back on pressure.', RuntimeWarning - ) + ) self.properties.ignition_target = 'pressure' self.properties.ignition_type = 'd/dt max' else: @@ -445,23 +406,27 @@ def run_case(self, restart=False): return # Save simulation results in hdf5 table format. - table_def = {'time': tables.Float64Col(pos=0), - 'temperature': tables.Float64Col(pos=1), - 'pressure': tables.Float64Col(pos=2), - 'volume': tables.Float64Col(pos=3), - 'mass_fractions': tables.Float64Col( - shape=(self.reac.thermo.n_species), pos=4 - ), - } - - with tables.open_file(self.meta['save-file'], mode='w', - title=self.meta['id'] - ) as h5file: - - table = h5file.create_table(where=h5file.root, - name='simulation', - description=table_def - ) + table_def = { + 'time': tables.Float64Col(pos=0), + 'temperature': tables.Float64Col(pos=1), + 'pressure': tables.Float64Col(pos=2), + 'volume': tables.Float64Col(pos=3), + 'mass_fractions': tables.Float64Col( + shape=(self.reac.thermo.n_species), pos=4 + ), + } + + with tables.open_file( + self.meta['save-file'], mode='w', + title=self.meta['id'] + ) as h5file: + + table = h5file.create_table( + where=h5file.root, + name='simulation', + description=table_def + ) + # Row instance to save timestep information to timestep = table.row # Save initial conditions @@ -491,7 +456,90 @@ def run_case(self, restart=False): # Write ``table`` to disk table.flush() - print('Done with case ', self.meta['id']) + print(f'Done with case {self.meta["id"]}') + + def get_ignition_delay(self, time, target): + """Function to get ignition delay + """ + ign_delays = [0.0] + # Analysis for ignition depends on type specified + if self.properties.ignition_type in ['max', 'd/dt max']: + if self.properties.ignition_type == 'd/dt max': + # Evaluate derivative + target = first_derivative(time.magnitude, target) + + # Get indices of peaks + ind = detect_peaks(target) + + # Fall back on derivative if max value doesn't work. + if len(ind) == 0 and self.properties.ignition_type == 'max': + target = first_derivative(time.magnitude, target) + ind = detect_peaks(target) + + # something has gone wrong if there is still no peak + if len(ind) == 0: + filename = 'target-data-' + self.meta['id'] + '.out' + warnings.warn( + 'No peak found, dumping target data to ' + + filename + ' and continuing', + RuntimeWarning + ) + np.savetxt( + filename, np.c_[time.magnitude, target], + header=('time, target (' + self.properties.ignition_target + ')') + ) + self.meta['simulated-ignition-delay'] = 0.0 * units.second + return + + # Get index of largest peak (overall ignition delay) + max_ind = ind[np.argmax(target[ind])] + + # Will need to subtract compression time for RCM + time_comp = 0.0 + if hasattr(self.properties.rcm_data, 'compression_time'): + if hasattr(self.properties.rcm_data.compression_time, 'value'): + time_comp = self.properties.rcm_data.compression_time.value + else: + time_comp = self.properties.rcm_data.compression_time + + ign_delays = time[ind[np.where( + (time[ind[ind <= max_ind]] - time_comp) + > 0. * units.second + )]] - time_comp + + elif self.properties.ignition_type == '1/2 max': + # maximum value, and associated index + max_val = np.max(target) + ind = detect_peaks(target) + max_ind = ind[np.argmax(target[ind])] + + # TODO: interpolate for actual half-max value + # Find index associated with the 1/2 max value, but only consider + # points before the peak + half_idx = (np.abs(target[0:max_ind] - 0.5 * max_val)).argmin() + ign_delays = time[half_idx] + + elif self.properties.ignition_type == 'd/dt max extrapolated': + # First need to evaluate derivative of the target + target_derivative = first_derivative(time, target) + + # Get indices of peaks, and index of largest peak, which corresponds to + # the point of maximum deriative + peak_inds = detect_peaks(target_derivative, edge=None, mph=1.e-9 * np.max(target)) + max_ind = peak_inds[np.argmax(target_derivative[peak_inds])] + + # use slope to extrapolate to intercept with baseline value (0 by default) + ign_delays = np.array([time.magnitude[max_ind] - (target[max_ind] / target_derivative[max_ind])]) * units.second + + else: + warnings.warn( + 'Unable to process ignition type ' + self.properties.ignition_type + + ', setting result to 0 and continuing', RuntimeWarning + ) + ign_delays = 0.0 * units.second + return np.array([0.0]) + + return ign_delays def process_results(self): """Process integration results to obtain ignition delay. @@ -509,34 +557,190 @@ def process_results(self): target = table.col('temperature') else: target = table.col('mass_fractions')[:, self.properties.ignition_target] - # store initial and maximum temperatures as a gate for whether the case ignited - init_temperature = table.col('temperature')[0] - max_temperature = np.max(table.col('temperature')) # add units to time time = time * units.second + ign_delays = self.get_ignition_delay(time, target) - # Will need to subtract compression time for RCM - time_comp = 0.0 - if hasattr(self.properties.rcm_data, 'compression_time'): - if hasattr(self.properties.rcm_data.compression_time, 'value'): - time_comp = self.properties.rcm_data.compression_time.value - else: - time_comp = self.properties.rcm_data.compression_time - - # First use basic check for ignition based on temperature increase of at least 50 K - if max_temperature >= init_temperature + 50.: - ignition_delays = get_ignition_delay(time.magnitude, target, - self.properties.ignition_target, - self.properties.ignition_type - ) - self.meta['simulated-ignition-delay'] = (ignition_delays[0] - time_comp) * units.second + # Overall ignition delay + if len(ign_delays) > 0: + self.meta['simulated-ignition-delay'] = ign_delays[-1] else: - warnings.warn('No ignition for case ' + self.meta['id'] + - ', setting value to 0.0 and continuing', - RuntimeWarning - ) self.meta['simulated-ignition-delay'] = 0.0 * units.second - # TODO: detect two-stage ignition. - self.meta['simulated-first-stage-delay'] = np.nan * units.second + # First-stage ignition delay + if len(ign_delays) > 1: + self.meta['simulated-first-stage-delay'] = ign_delays[0] + else: + self.meta['simulated-first-stage-delay'] = np.nan * units.second + + +class JSRSimulation(Simulation): + """Class for jet-stirred reaction simulation with the goal of + obtaining species concentration profiles""" + def __init__(self, *args, **kwargs): + if 'species_name' in kwargs: + self.target_species_name = kwargs['species_name'] + super(self.__class__, self).__init__(*args, **kwargs) + + def setup_case(self, model_file, species_key, path=''): + """Sets up the simulation case to be run. + + :param str model_file: Filename for Cantera-format model + :param dict species_key: Dictionary with species names for `model_file` + :param str path: Path for data file + """ + # Establishes the model + self.gas = ct.Solution(model_file) + self.species_key = species_key + # Set max simulation time, pressure valve coefficient, and max pressure rise for Cantera + # These could be set to something in ChemKED file, but haven't seen these specified at all.... + self.maxsimulationtime = 60 + self.pressurevalcof = 0.01 + self.maxpressurerise = 0.01 + # Reactor volume needed in m^3 for Cantera + self.volume = self.properties.reactor_volume.magnitude + # Residence time needed in s for Cantera + self.restime = self.properties.residence_time.magnitude + # Create reactants from chemked file + reactants = [ + self.species_key[self.properties.inlet_composition[spec].species_name] + ':' + + str(self.properties.inlet_composition[spec].amount.magnitude.nominal_value) + for spec in self.properties.inlet_composition + ] + + reactants = ','.join(reactants) + self.properties.pressure.ito('pascal') + + # Need to extract values from quantity or measurement object + if hasattr(self.properties.pressure, 'value'): + pres = self.properties.pressure.value.magnitude + elif hasattr(self.properties.pressure, 'nominal_value'): + pres = self.properties.pressure.nominal_value + else: + pres = self.properties.pressure.magnitude + temperature = self.properties.temperature + temperature.ito('kelvin') + if hasattr(temperature, 'value'): + temp = temperature.value.magnitude + elif hasattr(temperature, 'nominal_value'): + temp = temperature.nominal_value + else: + temp = temperature.magnitude + + if self.properties.inlet_composition_type in ['mole fraction', 'mole percent']: + self.gas.TPX = temp, pres, reactants + elif self.properties.inlet_composition_type == 'mass fraction': + self.gas.TPY = temp, pres, reactants + else: + raise(BaseException('error: not supported')) + return + # Upstream and exhaust + self.fuelairmix = ct.Reservoir(self.gas) + self.exhaust = ct.Reservoir(self.gas) + + # Ideal gas reactor + self.reactor = ct.IdealGasReactor(self.gas, energy='off', volume=self.volume) + self.massflowcontrol = ct.MassFlowController( + upstream=self.fuelairmix, + downstream=self.reactor, + mdot=self.reactor.mass / self.restime + ) + self.pressureregulator = ct.Valve( + upstream=self.reactor, + downstream=self.exhaust, + K=self.pressurevalcof + ) + + # Create reactor newtork + self.reactor_net = ct.ReactorNet([self.reactor]) + file_path = os.path.join(path, self.meta['id'] + '.h5') + self.meta['save-file'] = file_path + if self.target_species_name is None: + self.meta['target-species-index'] = -1 + else: + self.meta['target-species-index'] = self.gas.species_index(species_key[self.target_species_name]) + + def run_case(self, restart=False): + """Run simulation case set up ``setup_case``. + + :param bool restart: If ``True``, skip if results file exists. + """ + + if restart and os.path.isfile(self.meta['save-file']): + print('Skipped existing case ', self.meta['id']) + return + + # Save simulation results in hdf5 table format + table_def = { + 'time': tables.Float64Col(pos=0), + 'temperature': tables.Float64Col(pos=1), + 'pressure': tables.Float64Col(pos=2), + 'volume': tables.Float64Col(pos=3), + 'mole_fractions': tables.Float64Col(shape=(self.reactor.thermo.n_species), pos=4), + } + + with tables.open_file( + self.meta['save-file'], mode='w', + title=self.meta['id'] + ) as h5file: + + table = h5file.create_table( + where=h5file.root, + name='simulation', + description=table_def + ) + + # Row instance to save timestep information to + timestep = table.row + # Save initial conditions + timestep['time'] = self.reactor_net.time + timestep['temperature'] = self.reactor.T + timestep['pressure'] = self.reactor.thermo.P + timestep['volume'] = self.reactor.volume + timestep['mole_fractions'] = self.reactor.thermo.X + # Add ``timestep`` to table + timestep.append() + + integration_failed = False + # Main time integration loop; continue integration while time of + # the ``ReactorNet`` is less than specified end time. + while self.reactor_net.time < self.maxsimulationtime: + self.reactor_net.step() + + # Save new timestep information + timestep['time'] = self.reactor_net.time + timestep['temperature'] = self.reactor.T + timestep['pressure'] = self.reactor.thermo.P + timestep['volume'] = self.reactor.volume + timestep['mole_fractions'] = self.reactor.thermo.X + + # Add ``timestep`` to table + timestep.append() + if len(table) > 5000: + integration_failed = True + break + + if integration_failed: + # Failure is due to this bug + # https://github.com/Cantera/cantera/issues/1117 + warnings.warn( + 'Normal integration timed out in JSR simulation ' + self.meta['id'], + RuntimeWarning + ) + + # Write ``table`` to disk + table.flush() + + print('Done with case ', self.meta['id']) + + def process_results(self): + """ + Process integration results to obtain concentrations. + """ + # Load saved integration results + with tables.open_file(self.meta['save-file'], 'r') as h5file: + # Load Table with Group name simulation + table = h5file.root.simulation + concentration = table.col('mole_fractions')[:, self.meta['target-species-index']] + return concentration[-1] diff --git a/pyteck/tests/dataset_file_jsr.txt b/pyteck/tests/dataset_file_jsr.txt new file mode 100644 index 0000000..883f951 --- /dev/null +++ b/pyteck/tests/dataset_file_jsr.txt @@ -0,0 +1 @@ +testfile_jsr.yaml diff --git a/pyteck/tests/jsr_zhang_2016_phi_2.csv b/pyteck/tests/jsr_zhang_2016_phi_2.csv new file mode 100644 index 0000000..79ca666 --- /dev/null +++ b/pyteck/tests/jsr_zhang_2016_phi_2.csv @@ -0,0 +1,26 @@ +Temperature,n-heptane,oxygen,carbon monoxide,carbon dioxide,methane,acetylene,ethylene,ethane,propene,propane,propadiene,propyne,formaldehyde,ethylene oxide,acetaldehyde,1-butene,2-butene,"1,3-butadiene",n-butane,ethanol,propylene oxide,furan,acrolein,propanal,acetone,1-pentene,2-pentene,cyclopentene,"1,3-pentadiene",2-propen-1-ol,propan-1-ol,butenone,methyl-furan,butanal,2-butanone,1-hexene,2-butenal,acetic acid,pentanal + pentanone,benzene,dihydrofuran,3-heptene,2-heptene,1-heptene,2-hexanone,propanoic acid,2-ethyl-5-methyl-furan,2-methyl-dihydrofuranone,sum of cyclic ether with 5 atoms,sum of cyclic ether with 4 atoms,sum of cyclic ether with 3 atoms,2methanol-5methyl-tetrahydrofuranone,heptanone +500,5.07E-03,2.93E-02,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,1.06E-08,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,4.83E-07,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,4.91E-08,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00 +525,4.92E-03,2.86E-02,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,4.58E-06,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,7.63E-06,1.06E-08,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,3.81E-07,1.93E-06,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,4.91E-08,0.00E+00,2.51E-07,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,9.15E-07,0.00E+00,0.00E+00,0.00E+00,0.00E+00 +550,4.66E-03,2.85E-02,0.00E+00,0.00E+00,0.00E+00,0.00E+00,6.02E-06,0.00E+00,3.10E-06,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,3.66E-05,8.07E-07,0.00E+00,0.00E+00,0.00E+00,1.04E-06,7.55E-07,0.00E+00,0.00E+00,7.99E-05,6.61E-06,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,4.99E-07,1.07E-06,5.79E-06,1.11E-05,0.00E+00,0.00E+00,2.22E-05,9.05E-06,4.91E-08,3.49E-06,3.53E-06,1.48E-06,2.51E-07,4.16E-06,6.43E-06,2.41E-06,2.90E-06,4.60E-05,0.00E+00,2.85E-06,3.71E-07,1.01E-06 +575,4.16E-03,2.63E-02,2.43E-04,1.01E-04,1.46E-06,0.00E+00,3.08E-05,0.00E+00,1.70E-05,0.00E+00,0.00E+00,0.00E+00,6.01E-05,5.91E-06,2.30E-04,6.46E-06,0.00E+00,0.00E+00,0.00E+00,6.92E-06,3.22E-06,9.89E-07,9.00E-06,4.32E-04,3.05E-05,2.92E-06,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,7.19E-06,2.63E-06,3.18E-05,3.27E-05,1.86E-07,3.89E-06,7.87E-05,3.08E-05,4.91E-08,8.47E-06,2.03E-05,9.46E-06,6.20E-06,1.12E-05,4.15E-05,7.24E-06,3.52E-06,3.13E-04,1.00E-04,2.16E-05,4.30E-06,1.33E-05 +600,3.55E-03,2.33E-02,9.68E-04,2.51E-04,2.19E-06,0.00E+00,7.30E-05,6.92E-07,4.44E-05,0.00E+00,0.00E+00,0.00E+00,4.57E-04,1.06E-05,3.59E-04,1.78E-05,0.00E+00,4.21E-07,3.62E-07,1.04E-05,4.08E-06,8.40E-07,2.52E-05,6.65E-04,2.74E-05,6.50E-06,0.00E+00,0.00E+00,0.00E+00,1.07E-06,0.00E+00,1.73E-05,3.81E-06,5.65E-05,2.31E-05,7.91E-07,8.89E-06,5.47E-05,2.50E-05,4.91E-08,5.23E-06,3.66E-05,1.66E-05,9.57E-06,1.02E-05,1.78E-05,4.83E-06,0.00E+00,4.91E-04,2.06E-04,4.25E-05,8.05E-06,2.07E-05 +625,3.36E-03,2.31E-02,1.42E-03,2.67E-04,5.83E-06,0.00E+00,1.65E-04,1.38E-06,9.29E-05,0.00E+00,0.00E+00,0.00E+00,6.01E-04,1.30E-05,3.80E-04,3.55E-05,4.84E-07,1.05E-06,7.96E-07,8.18E-06,6.04E-06,2.62E-06,4.59E-05,7.74E-04,2.29E-05,1.48E-05,4.73E-07,4.73E-07,5.73E-07,1.83E-06,6.03E-07,3.04E-05,5.98E-06,6.76E-05,1.56E-05,1.86E-06,1.45E-05,1.82E-05,1.66E-05,4.91E-08,2.89E-06,5.58E-05,2.59E-05,1.04E-05,5.20E-06,1.29E-05,3.38E-06,0.00E+00,6.37E-04,3.06E-04,5.68E-05,9.17E-06,2.34E-05 +650,3.67E-03,2.45E-02,9.16E-04,1.46E-04,7.29E-06,0.00E+00,1.97E-04,1.04E-06,1.18E-04,0.00E+00,0.00E+00,0.00E+00,2.41E-04,5.91E-06,3.26E-04,5.21E-05,6.46E-07,1.89E-06,8.32E-07,2.75E-06,4.33E-06,2.08E-06,4.41E-05,6.76E-04,1.58E-05,2.57E-05,7.09E-07,1.03E-06,1.27E-06,2.33E-06,7.84E-07,1.65E-05,6.29E-06,6.51E-05,1.11E-05,2.84E-06,1.42E-05,6.20E-06,1.16E-05,4.91E-08,8.47E-07,7.10E-05,3.49E-05,1.02E-05,1.73E-06,7.91E-06,0.00E+00,0.00E+00,5.38E-04,3.26E-04,5.40E-05,8.08E-06,2.41E-05 +675,4.38E-03,2.77E-02,2.25E-04,0.00E+00,4.37E-06,0.00E+00,9.10E-05,0.00E+00,6.97E-05,0.00E+00,0.00E+00,0.00E+00,3.61E-05,2.36E-06,1.26E-04,4.08E-05,3.23E-07,1.26E-06,4.70E-07,1.64E-06,1.51E-06,6.43E-07,1.40E-05,3.01E-04,2.55E-06,2.66E-05,3.55E-07,5.61E-07,1.15E-06,8.18E-07,0.00E+00,3.99E-06,3.58E-06,3.47E-05,5.53E-06,3.56E-06,4.84E-06,0.00E+00,4.70E-06,4.91E-08,0.00E+00,5.65E-05,2.91E-05,8.36E-06,0.00E+00,0.00E+00,0.00E+00,0.00E+00,3.44E-04,2.36E-04,2.92E-05,5.06E-06,2.00E-05 +700,4.79E-03,2.87E-02,0.00E+00,0.00E+00,0.00E+00,0.00E+00,6.02E-06,0.00E+00,5.16E-06,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,1.26E-05,3.92E-06,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,2.29E-05,0.00E+00,3.37E-06,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,5.71E-07,3.47E-06,0.00E+00,5.58E-07,0.00E+00,0.00E+00,1.45E-06,4.91E-08,0.00E+00,4.78E-06,2.78E-06,7.52E-07,0.00E+00,0.00E+00,0.00E+00,0.00E+00,1.74E-05,9.88E-06,1.57E-06,3.61E-07,2.41E-06 +725,4.89E-03,2.93E-02,0.00E+00,0.00E+00,0.00E+00,0.00E+00,3.01E-06,0.00E+00,2.58E-06,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,5.72E-06,1.94E-06,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,1.05E-05,0.00E+00,1.68E-06,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,7.24E-07,0.00E+00,2.79E-07,0.00E+00,0.00E+00,0.00E+00,7.32E-08,0.00E+00,2.03E-06,1.15E-06,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,4.24E-06,0.00E+00,2.48E-07,0.00E+00,0.00E+00 +750,4.91E-03,2.84E-02,0.00E+00,0.00E+00,1.46E-06,0.00E+00,3.76E-06,0.00E+00,3.10E-06,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,3.43E-06,2.66E-06,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,9.08E-06,0.00E+00,2.25E-06,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,8.69E-07,0.00E+00,3.72E-07,0.00E+00,0.00E+00,0.00E+00,7.32E-08,0.00E+00,1.84E-06,9.61E-07,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,3.11E-06,0.00E+00,0.00E+00,0.00E+00,0.00E+00 +775,4.93E-03,2.80E-02,0.00E+00,0.00E+00,4.37E-06,0.00E+00,1.50E-05,0.00E+00,1.08E-05,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,6.86E-06,8.07E-06,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,1.27E-05,0.00E+00,5.32E-06,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,1.69E-06,0.00E+00,1.33E-06,0.00E+00,0.00E+00,0.00E+00,7.32E-08,0.00E+00,2.88E-06,1.78E-06,2.51E-07,0.00E+00,0.00E+00,0.00E+00,0.00E+00,5.06E-06,0.00E+00,0.00E+00,0.00E+00,0.00E+00 +800,4.78E-03,2.82E-02,0.00E+00,0.00E+00,1.46E-05,0.00E+00,5.19E-05,0.00E+00,2.84E-05,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,9.15E-06,2.30E-05,0.00E+00,4.21E-07,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,1.82E-05,0.00E+00,1.39E-05,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,7.62E-07,2.56E-06,0.00E+00,4.32E-06,0.00E+00,0.00E+00,0.00E+00,9.73E-08,0.00E+00,5.26E-06,2.95E-06,5.85E-07,0.00E+00,0.00E+00,0.00E+00,0.00E+00,7.56E-06,0.00E+00,0.00E+00,0.00E+00,0.00E+00 +825,4.41E-03,2.80E-02,1.49E-05,0.00E+00,4.37E-05,0.00E+00,1.91E-04,0.00E+00,8.78E-05,0.00E+00,0.00E+00,0.00E+00,0.00E+00,1.18E-06,1.60E-05,6.82E-05,0.00E+00,2.10E-06,0.00E+00,0.00E+00,0.00E+00,0.00E+00,1.44E-06,2.91E-05,0.00E+00,3.78E-05,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,1.60E-06,3.67E-06,0.00E+00,1.42E-05,0.00E+00,0.00E+00,2.17E-06,1.21E-07,0.00E+00,8.40E-06,5.72E-06,1.61E-06,0.00E+00,0.00E+00,0.00E+00,0.00E+00,1.20E-05,0.00E+00,8.26E-07,0.00E+00,0.00E+00 +850,3.68E-03,2.80E-02,4.18E-04,1.66E-04,2.60E-04,3.69E-07,1.26E-03,2.08E-05,4.91E-04,1.84E-06,5.38E-07,0.00E+00,8.90E-05,1.30E-05,7.32E-05,3.30E-04,1.57E-06,3.11E-05,2.32E-06,0.00E+00,3.22E-06,6.43E-07,2.07E-05,9.81E-05,0.00E+00,1.57E-04,2.07E-06,2.48E-06,5.73E-06,2.77E-06,0.00E+00,1.05E-06,5.33E-06,7.72E-06,0.00E+00,6.37E-05,3.49E-06,0.00E+00,3.62E-06,2.90E-07,0.00E+00,2.01E-05,1.21E-05,3.03E-06,0.00E+00,0.00E+00,0.00E+00,0.00E+00,2.46E-05,0.00E+00,2.48E-06,4.93E-07,3.13E-07 +875,2.13E-03,2.45E-02,1.65E-03,2.22E-04,6.33E-04,3.09E-06,3.05E-03,7.06E-05,9.67E-04,1.11E-05,1.34E-06,2.85E-06,3.61E-04,3.31E-05,1.58E-04,5.03E-04,4.04E-06,8.92E-05,5.79E-06,0.00E+00,7.55E-06,2.32E-06,9.00E-05,1.60E-04,9.15E-06,2.04E-04,6.79E-06,7.09E-06,1.12E-05,8.18E-06,7.84E-07,1.70E-06,4.95E-06,1.16E-05,4.52E-06,8.09E-05,1.13E-05,0.00E+00,4.70E-06,4.35E-07,0.00E+00,1.62E-05,9.53E-06,2.67E-06,0.00E+00,0.00E+00,0.00E+00,0.00E+00,1.45E-05,0.00E+00,1.36E-06,0.00E+00,0.00E+00 +900,1.03E-03,2.05E-02,5.51E-03,3.69E-04,1.24E-03,1.42E-05,5.23E-03,1.54E-04,1.19E-03,1.71E-05,2.55E-06,3.93E-06,5.77E-04,6.14E-05,2.69E-04,4.51E-04,5.65E-06,1.22E-04,7.96E-06,0.00E+00,8.56E-06,2.18E-06,1.46E-04,1.78E-04,2.03E-05,1.61E-04,9.16E-06,8.86E-06,1.12E-05,1.64E-05,1.15E-06,1.20E-06,2.74E-06,9.17E-06,9.05E-06,6.44E-05,1.75E-05,0.00E+00,3.62E-06,8.92E-07,0.00E+00,8.90E-06,4.80E-06,1.25E-06,0.00E+00,0.00E+00,0.00E+00,0.00E+00,5.90E-06,0.00E+00,6.61E-07,0.00E+00,0.00E+00 +925,5.82E-04,1.79E-02,8.59E-03,6.78E-04,1.62E-03,2.71E-05,5.78E-03,1.94E-04,1.00E-03,2.21E-05,2.89E-06,3.99E-06,6.01E-04,4.96E-05,2.61E-04,3.12E-04,5.25E-06,1.03E-04,7.60E-06,0.00E+00,5.54E-06,1.48E-06,1.19E-04,1.38E-04,1.78E-05,1.06E-04,8.27E-06,7.39E-06,7.54E-06,1.38E-05,1.27E-06,9.48E-07,1.14E-06,5.31E-06,7.04E-06,4.19E-05,1.45E-05,0.00E+00,2.17E-06,1.01E-06,0.00E+00,4.76E-06,2.59E-06,8.56E-07,0.00E+00,0.00E+00,0.00E+00,0.00E+00,2.43E-06,0.00E+00,0.00E+00,0.00E+00,0.00E+00 +950,3.88E-04,1.47E-02,1.05E-02,1.07E-03,1.89E-03,3.77E-05,5.74E-03,2.08E-04,7.81E-04,2.26E-05,3.09E-06,3.99E-06,5.77E-04,3.19E-05,2.15E-04,2.10E-04,4.44E-06,8.08E-05,6.51E-06,0.00E+00,2.47E-06,8.90E-07,8.73E-05,1.13E-04,1.27E-05,6.79E-05,6.79E-06,5.61E-06,5.13E-06,8.81E-06,9.05E-07,5.49E-07,5.33E-07,4.34E-06,4.52E-06,2.65E-05,8.64E-06,0.00E+00,1.81E-06,8.44E-07,0.00E+00,3.61E-06,1.48E-06,5.85E-07,0.00E+00,0.00E+00,0.00E+00,0.00E+00,1.18E-06,0.00E+00,0.00E+00,0.00E+00,0.00E+00 +975,2.35E-04,1.28E-02,1.19E-02,1.36E-03,1.98E-03,4.79E-05,5.12E-03,2.01E-04,5.53E-04,2.24E-05,2.69E-06,3.81E-06,3.37E-04,1.83E-05,1.57E-04,1.26E-04,3.43E-06,5.85E-05,4.89E-06,0.00E+00,7.55E-07,0.00E+00,4.37E-05,6.17E-05,3.57E-06,3.90E-05,4.43E-06,2.95E-06,3.02E-06,4.91E-06,6.03E-07,0.00E+00,0.00E+00,2.70E-06,2.41E-06,1.51E-05,4.79E-06,0.00E+00,1.45E-06,6.51E-07,0.00E+00,2.09E-06,6.06E-07,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00 +1000,1.14E-04,1.16E-02,1.34E-02,1.82E-03,2.11E-03,5.91E-05,4.45E-03,1.94E-04,3.41E-04,2.12E-05,2.02E-06,3.36E-06,1.44E-04,9.45E-06,1.01E-04,6.58E-05,2.62E-06,3.85E-05,3.80E-06,0.00E+00,0.00E+00,0.00E+00,2.21E-05,4.00E-05,0.00E+00,1.83E-05,2.75E-06,2.07E-06,1.81E-06,1.70E-06,0.00E+00,0.00E+00,0.00E+00,1.40E-06,0.00E+00,6.98E-06,2.00E-06,0.00E+00,1.09E-06,3.86E-07,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00 +1025,4.83E-05,9.88E-03,1.52E-02,2.41E-03,2.60E-03,6.89E-05,4.29E-03,1.94E-04,2.48E-04,1.98E-05,1.41E-06,2.73E-06,3.85E-05,6.14E-06,7.78E-05,3.55E-05,1.70E-06,2.90E-05,2.71E-06,0.00E+00,0.00E+00,0.00E+00,1.31E-05,3.09E-05,0.00E+00,7.98E-06,1.77E-06,1.30E-06,1.06E-06,6.29E-07,0.00E+00,0.00E+00,0.00E+00,9.65E-07,0.00E+00,3.02E-06,9.98E-07,0.00E+00,5.43E-07,1.70E-07,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00 +1050,1.64E-05,8.16E-03,1.83E-02,2.97E-03,3.15E-03,8.68E-05,4.11E-03,1.92E-04,1.49E-04,1.75E-05,9.41E-07,1.99E-06,1.54E-05,3.90E-06,4.58E-05,1.53E-05,1.01E-06,2.02E-05,1.59E-06,0.00E+00,0.00E+00,0.00E+00,4.05E-06,1.05E-05,0.00E+00,2.01E-06,7.39E-07,4.43E-07,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,4.83E-07,0.00E+00,6.98E-07,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00 +1075,1.22E-06,5.48E-03,1.95E-02,3.67E-03,2.46E-03,7.41E-05,2.31E-03,1.23E-04,4.34E-05,9.22E-06,0.00E+00,1.14E-06,0.00E+00,0.00E+00,1.83E-05,2.58E-06,4.04E-07,6.31E-06,5.07E-07,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,2.95E-07,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00 +1100,0.00E+00,3.24E-03,2.14E-02,4.38E-03,2.46E-03,8.24E-05,1.67E-03,9.76E-05,1.14E-05,4.61E-06,0.00E+00,0.00E+00,0.00E+00,0.00E+00,6.86E-06,8.07E-07,0.00E+00,2.10E-06,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00 diff --git a/pyteck/tests/nheptane_reduced.cti b/pyteck/tests/nheptane_reduced.cti new file mode 100644 index 0000000..36a68cb --- /dev/null +++ b/pyteck/tests/nheptane_reduced.cti @@ -0,0 +1,173 @@ +units(length='cm', time='s', quantity='mol', act_energy='kcal/mol') + +ideal_gas(name='gas', + elements="H D T C Ci O Oi N Ne Ar He Si S F Cl Br I X", + species="""N2 Ar He + C7H16(1) O2(2) + O(5) CO2(7) + H2(13) H(14) OH(15)""", + reactions='all', + transport='Mix', + initial_state=state(temperature=300.0, pressure=OneAtm)) + +#------------------------------------------------------------------------------- +# Element data +#------------------------------------------------------------------------------- + +element(symbol='Ci', atomic_mass=13.003) +element(symbol='D', atomic_mass=2.014) +element(symbol='Oi', atomic_mass=17.999) +element(symbol='T', atomic_mass=3.016) +element(symbol='X', atomic_mass=195.083) +#------------------------------------------------------------------------------- +# Species data +#------------------------------------------------------------------------------- + +species(name='N2', + atoms='N:2', + thermo=(NASA([100.00, 1817.05], + [ 3.61262543E+00, -1.00892648E-03, 2.49897620E-06, + -1.43374896E-09, 2.58633809E-13, -1.05110286E+03, + 2.65270528E+00]), + NASA([1817.05, 5000.00], + [ 2.97592894E+00, 1.64137380E-03, -7.19704298E-07, + 1.25374055E-10, -7.91499564E-15, -1.02585898E+03, + 5.53740798E+00])), + transport=gas_transport(geom='linear', + diam=3.621, + well_depth=97.53, + polar=1.76, + rot_relax=4.0)) + +species(name='Ar', + atoms='Ar:1', + thermo=(NASA([100.00, 4762.74], + [ 2.50000000E+00, -1.87769028E-11, 2.35957382E-14, + -9.62803665E-18, 1.20727680E-21, -7.45000000E+02, + 4.36630362E+00]), + NASA([4762.74, 5000.00], + [ 2.87626762E+00, -3.12580154E-04, 9.73654773E-08, + -1.34776059E-11, 6.99515521E-16, -1.10730235E+03, + 1.95965953E+00])), + transport=gas_transport(geom='atom', + diam=3.33, + well_depth=136.501)) + +species(name='He', + atoms='He:1', + thermo=(NASA([100.00, 4762.74], + [ 2.50000000E+00, -1.87769028E-11, 2.35957382E-14, + -9.62803665E-18, 1.20727680E-21, -7.45000000E+02, + 9.14221516E-01]), + NASA([4762.74, 5000.00], + [ 2.87626762E+00, -3.12580154E-04, 9.73654773E-08, + -1.34776059E-11, 6.99515521E-16, -1.10730235E+03, + -1.49242258E+00])), + transport=gas_transport(geom='atom', + diam=2.576, + well_depth=10.2)) + +species(name='C7H16(1)', + atoms='C:7 H:16', + thermo=(NASA([100.00, 1262.75], + [-2.94181528E-03, 7.69376326E-02, -3.48033705E-05, + 1.86184391E-09, 2.01542087E-12, -2.57829331E+04, + 3.00313064E+01]), + NASA([1262.75, 5000.00], + [ 1.38380723E+01, 4.85940646E-02, -1.95472187E-05, + 3.52831867E-09, -2.39060433E-13, -3.05142881E+04, + -4.48658921E+01])), + transport=gas_transport(geom='nonlinear', + diam=6.366, + well_depth=402.997)) + +species(name='O2(2)', + atoms='O:2', + thermo=(NASA([100.00, 1087.71], + [ 3.53763685E+00, -1.22827511E-03, 5.36758628E-06, + -4.93128119E-09, 1.45955119E-12, -1.03799025E+03, + 4.67179810E+00]), + NASA([1087.71, 5000.00], + [ 3.16427277E+00, 1.69453633E-03, -8.00335204E-07, + 1.59029939E-10, -1.14890928E-14, -1.04844563E+03, + 6.08302729E+00])), + transport=gas_transport(geom='linear', + diam=3.467, + well_depth=106.7)) + +species(name='O(5)', + atoms='O:1', + thermo=(NASA([100.00, 4762.74], + [ 2.50000000E+00, -1.87769028E-11, 2.35957382E-14, + -9.62803665E-18, 1.20727680E-21, 2.92267216E+04, + 5.11106769E+00]), + NASA([4762.74, 5000.00], + [ 2.87626762E+00, -3.12580154E-04, 9.73654773E-08, + -1.34776059E-11, 6.99515521E-16, 2.88644192E+04, + 2.70442360E+00])), + transport=gas_transport(geom='atom', + diam=2.75, + well_depth=80.0)) + +species(name='CO2(7)', + atoms='C:1 O:2', + thermo=(NASA([100.00, 988.19], + [ 3.27789791E+00, 2.75787976E-03, 7.12767589E-06, + -1.07852000E-08, 4.14216610E-12, -4.84756030E+04, + 5.97857463E+00]), + NASA([988.19, 5000.00], + [ 4.55073040E+00, 2.90725233E-03, -1.14641338E-06, + 2.25793556E-10, -1.69522790E-14, -4.89860167E+04, + -1.45672026E+00])), + transport=gas_transport(geom='linear', + diam=3.941, + well_depth=195.201)) + +species(name='H2(13)', + atoms='H:2', + thermo=(NASA([100.00, 1962.85], + [ 3.42253753E+00, 2.86648327E-04, -4.14667543E-07, + 4.27545515E-10, -9.38112761E-14, -1.02978491E+03, + -3.86364800E+00]), + NASA([1962.85, 5000.00], + [ 2.74217098E+00, 5.79561035E-04, 1.97192680E-07, + -6.41074329E-11, 4.95988421E-15, -5.52027985E+02, + 4.14195865E-01])), + transport=gas_transport(geom='linear', + diam=2.833, + well_depth=59.7)) + +species(name='H(14)', + atoms='H:1', + thermo=(NASA([100.00, 4762.74], + [ 2.50000000E+00, -1.87769028E-11, 2.35957382E-14, + -9.62803665E-18, 1.20727680E-21, 2.54727081E+04, + -4.59566260E-01]), + NASA([4762.74, 5000.00], + [ 2.87626762E+00, -3.12580154E-04, 9.73654773E-08, + -1.34776059E-11, 6.99515521E-16, 2.51104058E+04, + -2.86621035E+00])), + transport=gas_transport(geom='atom', + diam=2.05, + well_depth=145.0)) + +species(name='OH(15)', + atoms='H:1 O:1', + thermo=(NASA([100.00, 1005.26], + [ 3.48580083E+00, 1.33390756E-03, -4.70021880E-06, + 5.64351291E-09, -2.06305753E-12, 3.41195681E+03, + 1.99785890E+00]), + NASA([1005.26, 5000.00], + [ 2.88223155E+00, 1.03872344E-03, -2.35672354E-07, + 1.40277748E-11, 6.34181035E-16, 3.66956895E+03, + 5.59065082E+00])), + transport=gas_transport(geom='linear', + diam=2.75, + well_depth=80.0)) + +#------------------------------------------------------------------------------- +# Reaction data +#------------------------------------------------------------------------------- + +# Reaction 1 +reaction('O2(2) + H(14) <=> O(5) + OH(15)', [1.040000e+14, 0.0, 15.286]) diff --git a/pyteck/tests/species_keys_jsr.yaml b/pyteck/tests/species_keys_jsr.yaml new file mode 100644 index 0000000..78d86f6 --- /dev/null +++ b/pyteck/tests/species_keys_jsr.yaml @@ -0,0 +1,12 @@ +--- +nheptane_reduced.cti: + H2: "H2(13)" + oxygen: "O2(2)" + O2: "O2(2)" + Ar: "Ar" + nC7H16: "C7H16(1)" + n-heptane: "C7H16(1)" + N2: "N2" + nitrogen: "N2" + CO2: "CO2(7)" + He: "He" diff --git a/pyteck/tests/test_eval_model.py b/pyteck/tests/test_eval_model.py index f25f1a8..948c4da 100644 --- a/pyteck/tests/test_eval_model.py +++ b/pyteck/tests/test_eval_model.py @@ -9,7 +9,7 @@ # Third-party libraries import numpy import pytest -from pyked.chemked import ChemKED, DataPoint +from pyked.chemked import ChemKED, IgnitionDataPoint # Taken from http://stackoverflow.com/a/22726782/1569494 try: @@ -86,9 +86,10 @@ def test_normal_dist_noise(self): # add normally distributed noise, standard deviation of 1.0 noise = numpy.random.normal(0.0, 1.0, num) - standard_dev = eval_model.estimate_std_dev(changing_variable, - dependent_variable + noise - ) + standard_dev = eval_model.estimate_std_dev( + changing_variable, + dependent_variable + noise + ) assert numpy.isclose(1.0, standard_dev, rtol=1.e-2) def test_repeated_points(self): @@ -98,9 +99,10 @@ def test_repeated_points(self): dependent_variable = numpy.arange(1, 10) changing_variable[1] = changing_variable[0] - standard_dev = eval_model.estimate_std_dev(changing_variable, - dependent_variable - ) + standard_dev = eval_model.estimate_std_dev( + changing_variable, + dependent_variable + ) assert standard_dev == eval_model.min_deviation @@ -110,15 +112,15 @@ class TestGetChangingVariable: def test_single_point(self): """Check normal behavior for single point. """ - cases = [DataPoint({'pressure': [numpy.random.rand(1) * units('atm')], - 'temperature': [numpy.random.rand(1) * units('K')], - 'composition': - {'kind': 'mole fraction', - 'species': [{'species-name': 'O2', 'amount': [1.0]}] - }, - 'ignition-type': None - }) - ] + cases = [IgnitionDataPoint({ + 'pressure': [numpy.random.rand(1) * units('atm')], + 'temperature': [numpy.random.rand(1) * units('K')], + 'composition': { + 'kind': 'mole fraction', + 'species': [{'species-name': 'O2', 'amount': [1.0]}] + }, + 'ignition-type': None, + })] variable = eval_model.get_changing_variable(cases) assert len(variable) == 1 @@ -132,14 +134,15 @@ def test_temperature_changing(self): temperatures = numpy.random.rand(num) * units('K') cases = [] for temp in temperatures: - dp = DataPoint({'pressure': [str(pressure[0])], - 'temperature': [str(temp)], - 'composition': - {'kind': 'mole fraction', - 'species': [{'species-name': 'O2', 'amount': [1.0]}] - }, - 'ignition-type': None - }) + dp = IgnitionDataPoint({ + 'pressure': [str(pressure[0])], + 'temperature': [str(temp)], + 'composition': { + 'kind': 'mole fraction', + 'species': [{'species-name': 'O2', 'amount': [1.0]}] + }, + 'ignition-type': None + }) cases.append(dp) variable = eval_model.get_changing_variable(cases) @@ -155,14 +158,15 @@ def test_pressure_changing(self): temperature = numpy.random.rand(1) * units('K') cases = [] for pres in pressures: - dp = DataPoint({'pressure': [str(pres)], - 'temperature': [str(temperature[0])], - 'composition': - {'kind': 'mole fraction', - 'species': [{'species-name': 'O2', 'amount': [1.0]}] - }, - 'ignition-type': None - }) + dp = IgnitionDataPoint({ + 'pressure': [str(pres)], + 'temperature': [str(temperature[0])], + 'composition': { + 'kind': 'mole fraction', + 'species': [{'species-name': 'O2', 'amount': [1.0]}] + }, + 'ignition-type': None + }) cases.append(dp) variable = eval_model.get_changing_variable(cases) @@ -178,14 +182,15 @@ def test_both_changing(self): temperatures = numpy.random.rand(num) * units('K') cases = [] for pres, temp in zip(pressures, temperatures): - dp = DataPoint({'pressure': [str(pres)], - 'temperature': [str(temp)], - 'composition': - {'kind': 'mole fraction', - 'species': [{'species-name': 'O2', 'amount': [1.0]}] - }, - 'ignition-type': None - }) + dp = IgnitionDataPoint({ + 'pressure': [str(pres)], + 'temperature': [str(temp)], + 'composition': { + 'kind': 'mole fraction', + 'species': [{'species-name': 'O2', 'amount': [1.0]}] + }, + 'ignition-type': None + }) cases.append(dp) with pytest.warns(RuntimeWarning, @@ -196,6 +201,7 @@ def test_both_changing(self): assert len(variable) == num assert numpy.allclose(variable, [c.temperature.magnitude for c in cases]) + class TestEvalModel: """ """ @@ -218,7 +224,25 @@ def test(self): results_path=temp_dir, num_threads=1, skip_validation=True - ) - assert numpy.isclose(output['average error function'], 58.78211242028232, rtol=1.e-3) + ) + # average error was already failing before JSR addition + assert numpy.isclose(output['average error function'], 58.78211242028232, rtol=2.0e-3) assert numpy.isclose(output['error function standard deviation'], 0.0, rtol=1.e-3) assert numpy.isclose(output['average deviation function'], 7.635983785416241, rtol=1.e-3) + + def test_jsr(self): + """Test overall evaluation of JSR model. + """ + + with TemporaryDirectory() as temp_dir: + output = eval_model.evaluate_model( + model_name='nheptane_reduced.cti', + spec_keys_file=self.relative_location('species_keys_jsr.yaml'), + dataset_file=self.relative_location('dataset_file_jsr.txt'), + species_name='n-heptane', + data_path=self.relative_location(''), + model_path=self.relative_location(''), + results_path=temp_dir, + num_threads=2, + skip_validation=True + ) diff --git a/pyteck/tests/test_simulation.py b/pyteck/tests/test_simulation.py index 660bc29..e11043a 100644 --- a/pyteck/tests/test_simulation.py +++ b/pyteck/tests/test_simulation.py @@ -16,7 +16,7 @@ print("Error: Cantera must be installed.") raise -from pyked.chemked import ChemKED, DataPoint, TimeHistory +from pyked.chemked import ChemKED, TimeHistory, IgnitionDataPoint # Taken from http://stackoverflow.com/a/22726782/1569494 try: @@ -110,7 +110,7 @@ def test_sample_pressure_rise(self): assert times[-1] == time_end # Ensure final pressure correct, and check constant derivative - assert np.allclose(pressures[-1], pres*(pres_rise * time_end + 1)) + assert np.allclose(pressures[-1], pres * (pres_rise * time_end + 1)) dpdt = simulation.first_derivative(times, pressures) assert np.allclose(dpdt, pres * pres_rise) @@ -122,8 +122,8 @@ def test_volume_profile_no_pressure_rise(self): """Ensure constant volume history if zero pressure rise. """ [times, volume] = simulation.create_volume_history( - 'air.xml', 300., ct.one_atm, 'N2:1.0', 0.0, 1.0 - ) + 'air.xml', 300., ct.one_atm, 'N2:1.0', 0.0, 1.0 + ) # check that end time is correct and volume unchanged assert np.isclose(times[-1], 1.0) assert np.allclose(volume, 1.0) @@ -136,9 +136,9 @@ def test_artificial_volume_profile_nitrogen(self): end_time = 1.0 initial_temp = 300. [times, volumes] = simulation.create_volume_history( - 'air.xml', initial_temp, initial_pres, 'N2:1.0', - pres_rise, end_time - ) + 'air.xml', initial_temp, initial_pres, 'N2:1.0', + pres_rise, end_time + ) # pressure at end time end_pres = initial_pres * (pres_rise * end_time + 1.0) @@ -199,12 +199,12 @@ def test_artificial_volume_profile(self): velocity_profile = simulation.PressureRiseProfile( 'air.xml', init_temp, init_pressure, 'N2:1.0', pressure_rise, end_time - ) + ) # Sample pressure [times, pressures] = simulation.sample_rising_pressure( end_time, init_pressure, 2.e3, pressure_rise - ) + ) # Check velocity profile against "theoretical" volume derivative gas = ct.Solution('air.xml') @@ -216,8 +216,8 @@ def test_artificial_volume_profile(self): gas.SP = init_entropy, pressures[i] gamma = gas.cp / gas.cv velocities[i] = velocity_profile(times[i]) - dvolumes[i] = ((-1. / gamma) * pressure_rise * - (pressures[i] / init_pressure)**((-1. / gamma) - 1.0) + dvolumes[i] = ((-1. / gamma) * pressure_rise + * (pressures[i] / init_pressure)**((-1. / gamma) - 1.0) ) assert np.allclose(velocities, dvolumes, rtol=1e-3) @@ -234,12 +234,22 @@ def test_max_species(self): """ a, b, c = [5.13293528e+04, 3.16147043e-01, 1.05018205e-02] times = np.linspace(0, 1, 10000) - mass_fraction = a * np.exp(-((times - b)/c)**2) + mass_fraction = a * np.exp(-((times - b) / c)**2) # max value of this occurs when x == b - - ignition_delays = simulation.get_ignition_delay(times, mass_fraction, 'species', 'max') - - assert np.allclose(ignition_delays[0], b, rtol=1e-4) + temperature = 1000.0 + times *= units.second + properties = IgnitionDataPoint({ + 'target name': 'species', + 'temperature': [str(temperature)], + 'composition': { + 'kind': 'mole fraction', + 'species': [{'species-name': 'O2', 'amount': [1.0]}] + }, + 'ignition-type': 'max' + }) + sim = simulation.AutoIgnitionSimulation('ignition delay', 'shock tube', {}, properties) + ignition_delays = sim.get_ignition_delay(times, mass_fraction) + assert np.allclose(ignition_delays.magnitude[-1], b, rtol=1e-4) def test_max_derivative(self): """Test using maximum derivative of temperature for ignition delay. @@ -249,67 +259,127 @@ def test_max_derivative(self): times = np.linspace(0, 1, 10000) temperature = -0.5 * np.sqrt(np.pi) * a * c * erf((b - times) / c) + d # max derivative of this occurs when x == b - - ignition_delays = simulation.get_ignition_delay(times, temperature, 'temperature', 'd/dt max') - - assert np.allclose(ignition_delays[0], b, rtol=1e-4) + times *= units.second + properties = IgnitionDataPoint({ + 'temperature': [str(temperature)], + 'composition': { + 'kind': 'mole fraction', + 'species': [{'species-name': 'O2', 'amount': [1.0]}] + }, + 'ignition-type': 'd/dt max' + }) + sim = simulation.AutoIgnitionSimulation('ignition delay', 'shock tube', {}, properties) + ignition_delays = sim.get_ignition_delay(times, temperature) + assert np.allclose(ignition_delays.magnitude[-1], b, rtol=1e-4) def test_max_derivative_species(self): """Test using max derivative of a species-looking profile. """ a, b, c = [5.13293528e+04, 3.16147043e-01, 1.05018205e-02] times = np.linspace(0, 1, 10000) - mass_fraction = a * np.exp(-((times - b)/c)**2) + mass_fraction = a * np.exp(-((times - b) / c)**2) + temperature = 1000.0 # first inflection point of Gaussian occurs at b - sqrt(1/2)*c # so this is where the maximum derivative occurs - - ignition_delays = simulation.get_ignition_delay(times, mass_fraction, 'species', 'd/dt max') - assert np.allclose(ignition_delays[0], b - np.sqrt(1/2)*c, rtol=1e-4) - + times *= units.second + properties = IgnitionDataPoint({ + 'temperature': [str(temperature)], + 'composition': { + 'kind': 'mole fraction', + 'species': [{'species-name': 'O2', 'amount': [1.0]}] + }, + 'ignition-type': 'd/dt max', + 'target name': 'species', + }) + sim = simulation.AutoIgnitionSimulation('ignition delay', 'shock tube', {}, properties) + ignition_delays = sim.get_ignition_delay(times, mass_fraction) + # ignition_delays = simulation.get_ignition_delay(times, mass_fraction, 'species', 'd/dt max') + assert np.allclose(ignition_delays.magnitude[-1], b - np.sqrt(1 / 2) * c, rtol=1e-4) def test_half_max(self): """Test using half maximum value for ignition delay. """ a, b, c = [5.13293528e+04, 3.16147043e-01, 1.05018205e-02] times = np.linspace(0, 1, 10000) - mass_fraction = a * np.exp(-((times - b)/c)**2) + mass_fraction = a * np.exp(-((times - b) / c)**2) + temperature = 1000.0 # value of peak is `a`, so half max is `a/2` # `mass_fraction = a/2` at `b - c*np.sqrt(np.log(2))` # (peak minus half of the full width and half max, FWHM) - - ignition_delays = simulation.get_ignition_delay(times, mass_fraction, 'species', '1/2 max') - - assert np.allclose(ignition_delays[0], b - c*np.sqrt(np.log(2)), rtol=1e-4) + times *= units.second + properties = IgnitionDataPoint({ + 'temperature': [str(temperature)], + 'composition': { + 'kind': 'mole fraction', + 'species': [{'species-name': 'O2', 'amount': [1.0]}] + }, + 'ignition-type': '1/2 max', + 'target name': 'species', + }) + sim = simulation.AutoIgnitionSimulation('ignition delay', 'shock tube', {}, properties) + ignition_delays = sim.get_ignition_delay(times, mass_fraction) + assert np.allclose(ignition_delays.magnitude, b - c * np.sqrt(np.log(2)), rtol=1e-4) def test_derivative_max_extrapolated(self): """Test using d/dt max extrapolated value for ignition delay. """ a, b, c = [5.13293528e+04, 3.16147043e-01, 1.05018205e-02] times = np.linspace(0, 1, 10000) - mass_fraction = a * np.exp(-((times - b)/c)**2) + mass_fraction = a * np.exp(-((times - b) / c)**2) + temperature = 1000.0 # first inflection point of Gaussian occurs at b - sqrt(1/2)*c # so this is where the maximum derivative occurs # derivative: # df_dt = (-2*a/c**2) * (times - b) * np.exp(-(b - times)**2 / c**2) - time_max_dfdt = b - np.sqrt(1/2)*c - dfdt_max = (-2*a/c**2) * (time_max_dfdt - b) * np.exp(-(b - time_max_dfdt)**2 / c**2) - - ignition_delays = simulation.get_ignition_delay(times, mass_fraction, 'species', 'd/dt max extrapolated') - assert np.allclose(ignition_delays[0], - time_max_dfdt - a * np.exp(-((time_max_dfdt - b)/c)**2) / dfdt_max, - rtol=1e-4 - ) + time_max_dfdt = b - np.sqrt(1 / 2) * c + dfdt_max = (-2 * a / c**2) * (time_max_dfdt - b) * np.exp(-(b - time_max_dfdt)**2 / c**2) + + times *= units.second + properties = IgnitionDataPoint({ + 'temperature': [str(temperature)], + 'composition': { + 'kind': 'mole fraction', + 'species': [{'species-name': 'O2', 'amount': [1.0]}] + }, + 'ignition-type': 'd/dt max extrapolated', + 'target name': 'species', + }) + sim = simulation.AutoIgnitionSimulation('ignition delay', 'shock tube', {}, properties) + ignition_delays = sim.get_ignition_delay(times, mass_fraction) + + # ignition_delays = simulation.get_ignition_delay(times, mass_fraction, 'species', 'd/dt max extrapolated') + assert np.allclose( + ignition_delays.magnitude[-1], + time_max_dfdt - a * np.exp(-((time_max_dfdt - b) / c)**2) / dfdt_max, + rtol=1e-4 + ) def test_not_supported_type(self): """Test that a non-supported type raises a warning and returns zero. """ - with pytest.warns(RuntimeWarning, + times = np.linspace(0, 1, 10000) + mass_fraction = times + temperature = 1000.0 + times *= units.second + properties = IgnitionDataPoint({ + 'temperature': [str(temperature)], + 'composition': { + 'kind': 'mole fraction', + 'species': [{'species-name': 'O2', 'amount': [1.0]}] + }, + 'ignition-type': 'min', + 'target name': 'emission', + }) + sim = simulation.AutoIgnitionSimulation('ignition delay', 'shock tube', {}, properties) + with pytest.warns( + RuntimeWarning, match='Unable to process ignition type min, setting result to 0 and continuing' - ): - ignition_delays = simulation.get_ignition_delay(0.0, 0.0, 'emission', 'min') + ): + # ignition_delays = simulation.get_ignition_delay(0.0, 0.0, 'emission', 'min') + ignition_delays = sim.get_ignition_delay(times, mass_fraction) - assert ignition_delays[0] == 0.0 + assert ignition_delays[-1] == 0.0 class TestSimulation: @@ -440,10 +510,10 @@ def test_shock_tube_pressure_rise_setup_case(self): # Check constructed velocity profile [times, volumes] = simulation.create_volume_history( - mechanism_filename, init_temp, init_pres, - 'H2:0.00444,O2:0.00566,AR:0.9899', - 0.10 * 1000., sim.time_end - ) + mechanism_filename, init_temp, init_pres, + 'H2:0.00444,O2:0.00566,AR:0.9899', + 0.10 * 1000., sim.time_end + ) volumes = volumes / volumes[0] dVdt = simulation.first_derivative(times, volumes) velocities = np.zeros(times.size) @@ -517,7 +587,7 @@ def test_rcm_setup_case(self): 5.78058719368E+001, 5.80849611077E+001, 5.85928651155E+001, 5.94734357453E+001, 6.09310671165E+001, 6.32487551103E+001, 6.68100309742E+001 - ]) + ]) volumes = volumes / volumes[0] dVdt = simulation.first_derivative(times, volumes) velocities = np.zeros(times.size) @@ -581,7 +651,7 @@ def test_shock_tube_run_cases(self): 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 9.95297294e-01, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00 - ]) + ]) assert np.allclose(table.col('time')[-1], time_end) assert np.allclose(table.col('temperature')[-1], temp) assert np.allclose(table.col('pressure')[-1], pres) @@ -626,7 +696,7 @@ def test_shock_tube_run_cases(self): 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 9.95297294e-01, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00 - ]) + ]) assert np.allclose(table.col('time')[-1], time_end) assert np.allclose(table.col('temperature')[-1], temp) assert np.allclose(table.col('pressure')[-1], pres) @@ -689,7 +759,7 @@ def test_shock_tube_pressure_rise_run_cases(self): 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 9.95297294e-01, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00 - ]) + ]) assert np.allclose(table.col('time')[-1], time_end) assert np.allclose(table.col('temperature')[-1], temp) assert np.allclose(table.col('pressure')[-1], pres) @@ -725,9 +795,10 @@ def test_rcm_run_cases(self): table = h5file.root.simulation # Ensure exact columns present - assert set(['time', 'temperature', 'pressure', - 'volume', 'mass_fractions' - ]) == set(table.colnames) + assert set([ + 'time', 'temperature', 'pressure', + 'volume', 'mass_fractions' + ]) == set(table.colnames) # Ensure final state matches expected time_end = 1.0e-1 temp = 2385.3726323703772 @@ -751,7 +822,7 @@ def test_rcm_run_cases(self): -1.31976752e-30, -2.12060990e-32, 1.55792718e-01, 7.74803838e-01, 2.72630502e-66, 2.88273784e-67, -2.18774836e-50, -1.47465442e-48 - ]) + ]) assert np.allclose(table.col('time')[-1], time_end) assert np.allclose(table.col('temperature')[-1], temp, rtol=1e-5, atol=1e-9 diff --git a/pyteck/tests/testfile_jsr.yaml b/pyteck/tests/testfile_jsr.yaml new file mode 100644 index 0000000..bc9bb2c --- /dev/null +++ b/pyteck/tests/testfile_jsr.yaml @@ -0,0 +1,177 @@ +--- +file-authors: + - name: Anthony Stohr + ORCID: 0000-0002-2388-1723 + - name: Benjamin Hoare + ORCID: 0000-0001-7205-0777 +file-version: 1.0 +chemked-version: 0.4.1 +reference: + doi: 10.1016/j.combustflame.2016.06.028 + authors: + - name: Kuiwen Zhang + - name: Colin Banyon + - name: John Bugler + - name: Henry J. Curran + ORCID: 0000-0002-5124-8562 + - name: Anne Rodriguez + - name: Olivier Herbinet + - name: Frédérique Battin-Leclerc + ORCID: 0000-0001-8265-7492 + - name: Christine BChir + - name: Karl Alexander Heufer + ORCID: 0000-0003-1657-5231 + journal: Combustion and Flame + year: 2016 + volume: 172 + pages: 116-135 + detail: An updated experimental and kinetic modeling study of n-heptane oxidation + +experiment-type: species profile +apparatus: + kind: jet stirred reactor + institution: University of Lorraine, Nancy, France + facility: LRGP JSR + +datapoints: + - csvfile: jsr_zhang_2016_phi_2.csv + temperature: + - column-name: Temperature # values are provided in the CSV file + - units: K + - uncertainty-type: absolute + uncertainty: 3 K + pressure: + - 1.067 bar + - uncertainty-type: absolute + uncertainty: 0.01 bar # made up for testing + reactor-volume: + - 0.000095 m^3 + residence-time: + - 2 s + inlet-composition: + kind: mole fraction + species: + - species-name: n-heptane + InChI: 1S/C7H16/c1-3-5-7-6-4-2/h3-7H2,1-2H3 + amount: + - 0.005 + - uncertainty-type: relative + uncertainty: 0.05 + - species-name: He + InChI: 1S/He + amount: + - 0.9675 + - uncertainty-type: relative + uncertainty: 0.05 + - species-name: oxygen + InChI: 1S/O2/c1-2 + amount: + - 0.0275 + - uncertainty-type: relative + uncertainty: 0.05 + outlet-composition: + kind: mole fraction + species: + - species-name: n-heptane # species names must correspond with column headers in the CSV file + InChI: 1S/C7H16/c1-3-5-7-6-4-2/h3-7H2,1-2H3 + - species-name: oxygen + InChI: InChI=1S/O2/c1-2 + - species-name: carbon monoxide + InChI: 1S/CH2O/c1-2/h1H2 + - species-name: carbon dioxide + InChI: 1S/CO2/c2-1-3 + - species-name: methane + InChI: 1S/CH4/h1H4 + - species-name: acetylene + InChI: 1S/C2H2/c1-2/h1-2H + - species-name: ethylene + InChI: 1S/C2H4/c1-2/h1-2H2 + - species-name: ethane + InChI: 1S/C2H6/c1-2/h1-2H3 + - species-name: propene + InChI: 1S/C3H6/c1-3-2/h3H,1H2,2H3 + - species-name: propane + InChI: 1S/C3H8/c1-3-2/h3H2,1-2H3 + - species-name: propadiene + InChI: 1S/C3H4/c1-3-2/h1-2H2 + - species-name: propyne + InChI: 1S/C3H4/c1-3-2/h1H,2H3 + - species-name: formaldehyde + InChI: 1S/CH2O/c1-2/h1H2 + - species-name: ethylene oxide + InChI: 1S/C2H4O/c1-2-3-1/h1-2H2 + - species-name: acetaldehyde + InChI: 1S/C2H4O/c1-2-3/h2H,1H3 + - species-name: 1-butene + InChI: 1S/C4H8/c1-3-4-2/h3H,1,4H2,2H3 + - species-name: 2-butene + InChI: 1S/C4H8/c1-3-4-2/h3-4H,1-2H3 + - species-name: 1,3-butadiene + InChI: 1S/C4H6/c1-3-4-2/h3-4H,1-2H2 + - species-name: n-butane + InChI: 1S/C4H10/c1-3-4-2/h3-4H2,1-2H3 + - species-name: ethanol + InChI: 1S/C2H6O/c1-2-3/h3H,2H2,1H3 + - species-name: propylene oxide + InChI: 1S/C3H6O/c1-3-2-4-3/h3H,2H2,1H3 + - species-name: furan + InChI: 1S/C4H4O/c1-2-4-5-3-1/h1-4H + - species-name: acrolein + InChI: 1S/C3H4O/c1-2-3-4/h2-3H,1H2 + - species-name: propanal + InChI: 1S/C3H8O/c1-2-3-4/h4H,2-3H2,1H3 + - species-name: acetone + InChI: 1S/C3H6O/c1-3(2)4/h1-2H3 + - species-name: 1-pentene + InChI: 1S/C5H10/c1-3-5-4-2/h3H,1,4-5H2,2H3 + - species-name: 2-pentene + InChI: 1S/C5H10/c1-3-5-4-2/h3,5H,4H2,1-2H3 + - species-name: cyclopentene + InChI: 1S/C5H8/c1-2-4-5-3-1/h1-2H,3-5H2 + - species-name: 1,3-pentadiene + InChI: 1S/C5H8/c1-3-5-4-2/h3-5H,1H2,2H3 + - species-name: 2-propen-1-ol + InChI: 1S/C3H6O/c1-2-3-4/h2,4H,1,3H2 + - species-name: propan-1-ol + InChI: 1S/C3H8O/c1-2-3-4/h4H,2-3H2,1H3 + - species-name: butenone + InChI: 1S/C4H6O/c1-3-4(2)5/h3H,1H2,2H3 + - species-name: methyl-furan + InChI: 1S/C5H6O/c1-5-3-2-4-6-5/h2-4H,1H3 + - species-name: butanal + InChI: 1S/C5H6O/c1-5-3-2-4-6-5/h2-4H,1H3 + - species-name: 2-butanone + InChI: 1S/C4H8O/c1-3-4(2)5/h3H2,1-2H3 + - species-name: 1-hexene + InChI: 1S/C6H12/c1-3-5-6-4-2/h3H,1,4-6H2,2H3 + - species-name: 2-butenal + InChI: 1S/C4H6O/c1-2-3-4-5/h2-4H,1H3 + - species-name: acetic acid + InChI: 1S/C2H4O2/c1-2(3)4/h1H3,(H,3,4) + - species-name: pentanal + pentanone + InChI: 1S/C5H10O/c1-2-3-4-5-6/h5H,2-4H2,1H3 + - species-name: benzene + InChI: 1S/C6H6/c1-2-4-6-5-3-1/h1-6H + - species-name: dihydrofuran + InChI: 1S/C4H6O/c1-2-4-5-3-1/h1,3H,2,4H2 + - species-name: 3-heptene + InChI: 1S/C7H14/c1-3-5-7-6-4-2/h5,7H,3-4,6H2,1-2H3 + - species-name: 2-heptene + InChI: 1S/C7H14/c1-3-5-7-6-4-2/h3,5H,4,6-7H2,1-2H3 + - species-name: 1-heptene + InChI: 1S/C7H14/c1-3-5-7-6-4-2/h3H,1,4-7H2,2H3 + - species-name: 2-hexanone + InChI: 1S/C6H12O/c1-3-4-5-6(2)7/h3-5H2,1-2H3 + - species-name: propanoic acid + InChI: 1S/C3H6O2/c1-2-3(4)5/h2H2,1H3,(H,4,5) + - species-name: 2-ethyl-5-methyl-furan + InChI: 1S/C7H10O/c1-3-7-5-4-6(2)8-7/h4-5H,3H2,1-2H3 + - species-name: 2-methyl-dihydrofuranone + InChI: 1S/C5H8O2/c1-4-5(6)2-3-7-4/h4H,2-3H2,1H3 + - species-name: sum of cyclic ether with 5 atoms + - species-name: sum of cyclic ether with 4 atoms + - species-name: sum of cyclic ether with 3 atoms + - species-name: 2methanol-5methyl-tetrahydrofuranone + InChI: 1S/C6H10O3/c1-4-6(8)2-5(3-7)9-4/h4-5,7H,2-3H2,1H3 + - species-name: heptanone + InChI: 1S/C7H14O/c1-3-4-5-6-7(2)8/h3-6H2,1-2H3