Negative CVs from NASA7CurveFit

[rebeccaem]

I know some of you already know about this issue, but I thought I'd document it and consolidate the threads here.

I'm trying to run methane combustion using a reaction mechanism from GRIMech. This includes some species that are not in the default Antioch files. I augmented the default species file (name, molecular weight, enthalpy of formation, etc) and I pulled the thermo data from GRIMech. This GRIMech thermo file is in the NASA7 format. With help from @pbauman, I set this up using the NASAThermoMixture, NASAEvaluator, etc.

It compiles and runs up to a point-- that is, I start at T_0, and the temperature increases for a bit until some T' at which point it drops suddenly to <1. Obviously this should not happen. I tried printing out some things to see what was going wrong and I noticed that some of the species return negative cv's.

Even without the kinetics, I still get the negative values. Here is a code snippet that produces them (I left the CEA lines that I had been using before in case it is helpful):

  Antioch::ChemicalMixture<double> chem_mixture( species_str_list, true, "speciesHCO.dat");
  //Antioch::CEAThermoMixture<double> cea_mixture( chem_mixture );
  Antioch::NASAThermoMixture<double, Antioch::NASA7CurveFit<double> > nasa_mixture( chem_mixture );
  Antioch::ReactionSet    <double> reaction_set( chem_mixture );

  Antioch::read_reaction_set_data_chemkin<double>( input_name, true, reaction_set );
  //Antioch::read_cea_mixture_data_ascii( cea_mixture, Antioch::DefaultFilename::thermo_data() );
  Antioch::read_nasa_mixture_data( nasa_mixture, "gri_thermo.dat", Antioch::CHEMKIN);

  Antioch::NASAEvaluator<double, Antioch::NASA7CurveFit<double> > thermo(nasa_mixture);

  typedef Antioch::TempCache<double> Cache;

  std::vector<double> h_RT_minus_s_R(n_species);
  thermo.h_RT_minus_s_R(Cache(temp),h_RT_minus_s_R);

  std::vector<double> cv(n_species, 0.);
  for (unsigned int s = 0; s < n_species; s++){
    cv[s] = thermo.cv(Cache(temp),s)/chem_mixture.R(s)* Antioch::Constants::R_universal<double>();
    if (cv[s] < 0.)
      std::cout<<"T is "<<temp<< ", species is " << s << " and cv  = "<< cv[s]<<"\n";
  }

For example, when temp = 1857.0, the cv of H2 as calculated above is -80.1385.

[pbauman]

@beanSnippet thanks for reporting. I've used a modified version of your program (plus a matplotlib script) to generate the images attached (I tried scalable pdfs, but GitHub only allows png, gif or jpg... sigh).

Something is definitely funky. @beanSnippet could you please do one additional step for me: could you take one of the obviously wrong species, e.g. HCCOH, and manually plot the polynomials for the coefficients in the .dat file (googling chemkin thermo will give you the polynomial form and the order of the coefficients, I don't remember it off the top of my head). Just something quick and dirty in Octave or something.

[variscarey]

Could you attach the gri_thermo data file? What I suspect is that (for
some reason, and this should be fixed, as it works elsewhere) that if
everything is correct in the input, the nasafit is not switching to the
high temperature polynomial at T=1000. This would be a bug on our part.
We may not be parsing the high temperature polynomial correctly. This
popped up during transport testing but I thought this was fixed.

Varis

On Tue, Aug 25, 2015 at 9:13 AM, Rebecca Morrison notifications@github.com
wrote:

I know some of you already know about this issue, but I thought I'd
document it and consolidate the threads here.

I'm trying to run methane combustion using a reaction mechanism from
GRIMech. This includes some species that are not in the default Antioch
files. I augmented the default species file (name, molecular weight,
enthalpy of formation, etc) and I pulled the thermo data from GRIMech. This
GRIMech thermo file is in the NASA7 format. With help from @pbauman
https://github.com/pbauman, I set this up using the NASAThermoMixture,
NASAEvaluator, etc.

It compiles and runs up to a point-- that is, I start at T_0, and the
temperature increases for a bit until some T' at which point it drops
suddenly to <1. Obviously this should not happen. I tried printing out some
things to see what was going wrong and I noticed that some of the species
return negative cv's.

Even without the kinetics, I still get the negative values. Here is a code
snippet that produces them (I left the CEA lines that I had been using
before in case it is helpful):

Antioch::ChemicalMixture chem_mixture( species_str_list, true, "speciesHCO.dat");
//Antioch::CEAThermoMixture cea_mixture( chem_mixture );
Antioch::NASAThermoMixture<double, Antioch::NASA7CurveFit > nasa_mixture( chem_mixture );
Antioch::ReactionSet reaction_set( chem_mixture );

Antioch::read_reaction_set_data_chemkin( input_name, true, reaction_set );
//Antioch::read_cea_mixture_data_ascii( cea_mixture, Antioch::DefaultFilename::thermo_data() );
Antioch::read_nasa_mixture_data( nasa_mixture, "gri_thermo.dat", Antioch::CHEMKIN);

Antioch::NASAEvaluator<double, Antioch::NASA7CurveFit > thermo(nasa_mixture);

typedef Antioch::TempCache Cache;

std::vector h_RT_minus_s_R(n_species);
thermo.h_RT_minus_s_R(Cache(temp),h_RT_minus_s_R);

std::vector cv(n_species, 0.);
for (unsigned int s = 0; s < n_species; s++){
cv[s] = thermo.cv(Cache(temp),s)/chem_mixture.R(s)* Antioch::Constants::R_universal();
if (cv[s] < 0.)
std::cout<<"T is "<<temp<< ", species is " << s << " and cv = "<< cv[s]<<"\n";
}

For example, when temp = 1857.0, the cv of H2 as calculated above is
-80.1385.

—
Reply to this email directly or view it on GitHub
#165.

Varis Carey
Assistant Professor
Department of Mathematical and Statistical Sciences
University of Colorado Denver

[pbauman]

@variscarey I'm emailing it to you privately as GitHub only allows posting images...

[pbauman]

@variscarey I figure it's probably something along the lines you mentioned. I'm hoping that @beanSnippet can manually do one of the polynomial sets and we can visually confirm that something like that is going on and focus our efforts. (Hopefully I'll finish #164 today so it will give us an easy way to add tests in for this wherever the bug is.)

[variscarey]

It's not switching to the right curve. I used a gri_thermo database I had
for Exact and the lower curve returns negative cv/R for high temperature.
(The almost linear curve is the high temperature cv/R)

Varis

On Tue, Aug 25, 2015 at 9:35 AM, Paul T. Bauman notifications@github.com
wrote:

@variscarey https://github.com/variscarey I figure it's probably
something along the lines you mentioned. I'm hoping that @beanSnippet
https://github.com/beanSnippet can manually do one of the polynomial
sets and we can visually confirm that something like that is going on and
focus our efforts. (Hopefully I'll finish #164
#164 today so it will give us
an easy way to add tests in for this wherever the bug is.)

—
Reply to this email directly or view it on GitHub
#165 (comment).

Varis Carey
Assistant Professor
Department of Mathematical and Statistical Sciences
University of Colorado Denver

[pbauman]

Thanks @variscarey. I'll try and look into this tonight, but someone feel free to beat me to it.

[pbauman]

Just FYI, I quickly glanced at the code and it the jump based on the interval is there, so it might be related to the parsing of multiple intervals. (Also, what is up with looping over the interval there... we shouldn't need to do that since we already grabbed what interval we're in...)

[pbauman]

(Also, what is up with looping over the interval there... we shouldn't need to do that since we already grabbed what interval we're in...)

Derp, nevermind. It's for when T is a vector type.

[rebeccaem]

It seems like you guys found the problem, but here's a plot! This is C_p/R vs T for H2.

Would it help to see H and S too?

[rebeccaem]

The jump is at 1000K here.

[rebeccaem]

One more point: in the thermo file, the coefficients for T>1000K are listed first, followed by T<1000K. At first I didn't catch that and I got a plot that looked similar to @pbauman's C_p/R plot above.

[SylvainPlessis]

@beanSnippet: this is probably the problem.

In the chemkin parser (file src/parsing/src/chemkin_parser, line 959 to 964; there), the a coeffs should be taken from 7 to 14 then from 0 to 6.

Basically, changing

for(unsigned int i = 0; i < 14; i++)
          {
            NumericType a;
            tmp >> a;
            coeffs.push_back(a);
          }

to

for(unsigned int i = 0; i < 7; i++)
          {
            NumericType a;
            tmp >> a;
            coeffs_high.push_back(a);
          }
for(unsigned int i = 0; i < 7; i++)
          {
            NumericType a;
            tmp >> a;
            coeffs_low.push_back(a);
          }
for(unsigned int i = 0; i < 14; i++)
          {
            (i < 7)?coeffs.push_back(coeffs_low[i]):coeffs.push_back(coeffs_high[i-7]);
          }

something of that sort.

Can someone try it quickly and verify this is the problem?

[SylvainPlessis]

Ok I had time to check and see, I indeed switch from negative values to positive correct-looking values with this bugfix. The chemkin parser indeed really needs better testing...

Thanks @beanSnippet for your remark.

[rebeccaem]

Thanks for fixing and verifying!

[pbauman]

#166 should fix this so I'm closing. @beanSnippet thanks for the report and please let us know if you encounter further issues. Regardless, hopefully we'll get better coverage of ChemKin parsing soon.