search

ckwatson / kernel

Does the math.
GNU General Public License v3.0
0 stars 0 forks source link

Concentrations all become NaN - #1

Closed ngraymon closed 8 years ago

ngraymon commented 8 years ago

For some reason the concentrations of the bromine puzzle all go to NaN. Plotting fails because of this.

ngraymon commented 8 years ago

So using a test case of 5mol of H2 and 0mol of Br2 with only a single H2 <-> H + H reaction this is what I observed.

The following error is key: /Users/ngraymon/Documents/projects/ckwatson/web_gui/kernel/engine/experiment_class.py:134 RuntimeWarning: divide by zero encountered in true_divide scaling_factor = np.divide(self.scaling_factor, np.amax(delta_energy_array)) /Users/ngraymon/Documents/projects/ckwatson/web_gui/kernel/engine/experiment_class.py:136 RuntimeWarning: invalid value encountered in multiply self.species_energy_array = np.multiply(scaling_factor, self.species_energy_array)

The call to get_reaction_rate_array(conc)(on line 275) returns a reaction_array array([ nan]) and therefore dYdt (on line 276) is array([ nan, nan]) and of course all of the concentrations afterwards will become NaN.

One possible option here would be to preform numerical checks inside the condition_elementary(conc, time) function. Although this would likely have an impact on performance.

Also we can see that the NaN values originate from the division by zero on line 134 scaling_factor = np.divide(self.scaling_factor, np.amax(delta_energy_array))

Investigating further it seems that the delta_energy_array (calculated on line 132) was zero

We look towards the product_energy_array and reactant_energy_array , which both turn out to be zero as a result of the species_energy_array being NaN. This should not be the case!

The species_energy_array is obtained (on line 96) as follows input_reaction_mechanism.get_energy_set()

Therefore the puzzle is failing because the reaction_mechanism has incorrect Gibbs energies for its chemical species.

It seems this is due to the incorrect Gibbs energy for hydrogen. It was 0.0, but it should be 203300.0. Once I fixed that (ckwatson/puzzles@e854c43566bd13c175db71fe1250dab48ea38ba8) the puzzle runs fine and generates plots I've added a small worst case scenario check to the engine. 914bd73d73d96c93505ea4844be0db24faf21710