diff --git a/t3/utils/flux.py b/t3/utils/flux.py index 6ef34655..903614fb 100644 --- a/t3/utils/flux.py +++ b/t3/utils/flux.py @@ -304,13 +304,24 @@ def set_batch_p(gas: ct.Solution, network.rtol = r_tol return network, reactor -def closest_bigger_number(array, target): - # Filter the numbers greater than the target +def closest_bigger_number(array: List[float], + target: float, + ) -> Tuple[Optional[float], Optional[int]]: + """ + Find the min number in an array that is greater than the target number. + Args: + - array (List[float]): List of numbers to search through. + - target (float): The target number. + Returns: + - Tuple[num[float], i[int]]: A tuple containing the min bigger number (num) and its index (i) in the originaul array. + If no number greater than the target is found, both elements in the tuple will be None. + """ + # Filter out the numbers smaller than the target candidates = [(num, i) for i, num in enumerate(array) if num > target] - # Find the closest bigger number + # Find the min bigger number if candidates: - m, i = min(candidates, key=lambda x: x[0]) - return m, i + num, i = min(candidates, key=lambda x: x[0]) + return num, i else: return None, None # No number greater than the target was found @@ -322,7 +333,7 @@ def run_jsr(gas: ct.Solution, V: Optional[float] = None, a_tol: float = 1e-16, r_tol: float = 1e-10, - tau_tolerance: float = 0.1, + tau_tolerance: float = 0.005, ) -> Dict[float, dict]: """ Run a JSR reactor with constant pressure, constant temperature, and constant volume. @@ -336,7 +347,9 @@ def run_jsr(gas: ct.Solution, V (float, optional): The reactor volume in cm^3. a_tol (float, optional): The absolute tolerance for the simulation. r_tol (float, optional): The relative tolerance for the simulation. - tau_tolerance (float, optional): Cantera will solve JSR equations with step() till tau*tolerance time, after that, it will solve them by advance(t_i). This will adress stiff equestion in which small time steps are required. + tau_tolerance (float, optional): Cantera will solve JSR equations with step() + till tau * tolerance time, after that, it will solve them by advance(t_i). + This will adress stiff equestion in which small time steps are required. Returns: dict: The T, P, X, and ROP profiles (values) at a specific time. @@ -347,34 +360,29 @@ def run_jsr(gas: ct.Solution, stoichiometry = get_rxn_stoichiometry(gas) for tau in times: t = 0 - dt = tau/num_steps + dt = tau / num_steps time_steps_array = np.arange(t, tau+dt, dt) #the last number is tau network, reactor = set_jsr(gas=gas, time=tau, composition=composition, T=T, P=P, V=V, a_tol=a_tol, r_tol=r_tol) - while t <= tau*tau_tolerance: + while t <= tau * tau_tolerance: network.step() - #update the system's time - t = network.time - if t > tau*tau_tolerance: - #get the current time of the system after step() - #find the closest_bigger time value in time_array, and continue with advance func t = network.time - t_array, t_step_index = closest_bigger_number(time_steps_array, t) - t = t_array - while t_step_index < len(time_steps_array) -1 : #till t =tau - rops = {spc.name: dict() for spc in gas.species()} - network.advance(t) - t_step_index +=1 - t = time_steps_array[t_step_index] - cantera_reaction_rops = gas.net_rates_of_progress - for spc in gas.species(): - for i, rxn in enumerate(gas.reactions()): - if stoichiometry[spc.name][i]: - if rxn.equation not in rops[spc.name].keys(): - rops[spc.name][rxn.equation] = 0 - rops[spc.name][rxn.equation] += cantera_reaction_rops[i] * stoichiometry[spc.name][i] - profile = {'P': gas.P, 'T': gas.T, 'X': {s.name: x for s, x in zip(gas.species(), gas.X)}, 'ROPs': rops} - if t == tau: - profiles[t] = profile + t = network.time + t, t_step_index = closest_bigger_number(time_steps_array, t) + while t_step_index < len(time_steps_array) -1 : #till t =tau + rops = {spc.name: dict() for spc in gas.species()} + network.advance(t) + t_step_index +=1 + t = time_steps_array[t_step_index] + cantera_reaction_rops = gas.net_rates_of_progress + for spc in gas.species(): + for i, rxn in enumerate(gas.reactions()): + if stoichiometry[spc.name][i]: + if rxn.equation not in rops[spc.name].keys(): + rops[spc.name][rxn.equation] = 0 + rops[spc.name][rxn.equation] += cantera_reaction_rops[i] * stoichiometry[spc.name][i] + profile = {'P': gas.P, 'T': gas.T, 'X': {s.name: x for s, x in zip(gas.species(), gas.X)}, 'ROPs': rops} + if t == tau: + profiles[t] = profile return profiles def run_batch_p(gas: ct.Solution,