ReactionMechanismGenerator / RMG-Java

The Java version of RMG: Reaction Mechanism Generator
http://rmg.sourceforge.net/
MIT License
29 stars 36 forks source link

Aborts when exceeds temperature and pressure range of rate coefficient. #10

Closed rwest closed 14 years ago

rwest commented 14 years ago

When this happens, RMG just reports Attempted to evaluate a rate coefficient outside the allowed temperature and pressure range. and stops. This is not a good failure mode!

rwest commented 14 years ago

In order of preference:

  1. figure out how to carry on
  2. carry on with the reaction systems that are still valid (if running several T,P simulations)
  3. report which rate, and which T,P, so that you can at least change the condition file or database!
rwest commented 14 years ago

Comment from mike at http://github.com/GreenGroup/RMG-Java/issues/closed#issue/11/comment/151211

Does this message occur when you are running a pressure-dependent job? If so, I'm guessing the temperature you are attempting to model is outside of 300 < T < 2100 K and/or the pressure is outside of 0.01 < P < 100 bar.

The following is speculation, but I know evaluating a Chebyshev pressure-dependent rate expression involves evaluating an arccos. If the T / P are outside of the range mentioned above, the argument to the arccos will be > 1 or < -1, which is not in the physical domain of arccos.

In these instances, should we tell the user to specify a different T / P range (which they can already do in the condition file) for which to compute the Chebyshev polynomials?

rwest commented 14 years ago

Yes, it's a pressure dependent job, but none of the conditions requested are outside the range you specify. Perhaps one is outside the range of a primary or seed reaction though. Relevant bits of condition file are at http://gist.github.com/322458

mrharper commented 14 years ago

I think I see the problem: You are asking for a PressureModel at 10 atm ... however your PRange only goes up to 10 bar.

rwest commented 14 years ago

Thanks. Trying again with 10 bar instead of 10 atm.

I think such inconsistencies should

  1. be caught earlier,
  2. be reported in such a way that they're easier to fix.
rwest commented 14 years ago

Failed again with 10 bar trying 9.5

rwest commented 14 years ago

Doesn't help. I've extended the error message a little, and am getting

Tried to evaluate P-dep rate coefficient at T=800.0K but only valid from 806.1 to 1963.1K
Reaction: O2(7)+C2H5(30)(+m)=SPC(38)+O(9)(+m)   1.0E0 0.0 0.0
TCHEB / 806.0629039746278 1963.0860015312323 /  PCHEB / 1.1515178797534147 8.458553284417798 /
CHEB / 4    3 /
CHEB / -7.2303000e-01 2.4402000e-16 1.7981000e-16 /
CHEB / 5.6453000e+00 7.7675000e-16 -7.1180000e-17 /
CHEB / 3.6916000e-02 4.9311000e-17 4.2178000e-16 /
CHEB / 5.9239000e-03 3.5068000e-16 -4.2554000e-17 /

For some reason the chebyshev polynomials, and the array of conditions at which k has been evaluated, do not span the requested range (800-2000K, 1-10bar):

this.chebyshev.Tlow = 806.063
this.chebyshev.Tup  = 1963.09
this.chebyshev.Plow = 116678
this.chebyshev.Pup  = 857063

These correspond to the lowest and highest of the this.temperatures and this.pressures arrays.

Any thoughts?

Stack trace:

0   0x00000002 in jing.rxn.PDepRateConstant.calculateRate() at PDepRateConstant.java:169
1   0x0000000f in jing.rxn.PDepReaction.calculateRate() at PDepReaction.java:341
2   0x000000b5 in jing.rxnSys.JDAS.transferReaction() at JDAS.java:573
3   0x0000020e in jing.rxnSys.JDAS.generatePDepODEReactionList() at JDAS.java:195
4   0x00000130 in jing.rxnSys.JDASSL.solve() at JDASSL.java:161
5   0x00000080 in jing.rxnSys.ReactionSystem.solveReactionSystem() at ReactionSystem.java:1500
6   0x00000283 in jing.rxnSys.ReactionModelGenerator.modelGeneration() at ReactionModelGenerator.java:1504
7   0x00000018 in RMG.main() at RMG.java:58

Unfortunately this only occurs 2.5 hours into the run, so debugging is a pain.

jwallen commented 14 years ago

Good news is I can reproduce the bug in < 1 min using a slightly modified version of the 1,3-hexadiene example network. The bug occurs after the TRange and PRange lines are read from the condition file, when the code chooses the exact set of temperatures and pressures to use. See line 4300 or so in ReactionModelGenerator.java.

From what I can tell, tempValueTilda and pressValueTilda are supposed to range from 1 to -1, but when I stepped through the code they only ranged from about 0.97 to -0.97, hence the bug. I have a fix in mind but should probably double-check with Mike before committing.

We should also add comments about why we distribute temperatures and pressures in this way.

mrharper commented 14 years ago

If you have access to the Chemkin 4.1.1 manual, the following discussion refers to the CHEMKIN_Theory documentation (if you have an MIT license, you can click the link below), in particular, near the bottom of page 49.

For determining which temperatures / pressures to calculate the rate coefficient at, in order to fit to a Chebyshev polynomial format, the Chemkin manual recommends using temperatures / pressures on the Gauss-Chebyshev grid. This is "to ensure the approximation is uniform over the desired domain".

As Josh mentioned, this results in the "normalized" temperatures / pressures being sent to fame to range from -0.97 to 0.97 (instead of -1 to 1). The problem is that although Mike thought he was setting the Tmin to 800K (using Richard's example), the source code must have a check to use the lowest temperature in the "temperatures" array as the Tlow; thus the 800K is changed to ~806K.

Possible solutions (w/o having discussed these w/Josh): 1) Use linear interpolation instead of points on the Gauss-Chebyshev grid 2) Introduce new methods to set the Tlow/Tup to be the user-specified values (and continue using the Gauss-Chebyshev grid points)

https://web.mit.edu/afs/athena/software/chemkin_v4.1.1/distrib/docs/CHEMKIN_Theory.pdf

jwallen commented 14 years ago

I propose we instead use Ttilde = cos(iPI/(dT-1)) for 0 <= i < dT and Ptilde = cos(jPI/(dP-1)) for 0 <= j < dP. This would give the proper endpoints and, if I'm understanding the Chemkin manual correctly, still give a Gauss-Chebyshev grid.

mrharper commented 14 years ago

Good suggestion Josh. Seems valid/reasonable to me.

rwest commented 14 years ago

I just pushed a bunch of changes to master that help with error messages, and make this sort of thing easier to catch.

rwest commented 14 years ago

Hopefully closed by 0260e398517e53b766dce9e80738bd53833d8ce1 Change Gauss-Chebyshev Grid definition

Previously it didn't actually reach the endpoints, so if you asked for Temperatures from 800-2000K you may get a grid from something like 806.5 to 1963.5.

Hopefully these nodes span the full range and are still roots of the Chebyshev polynomials of the first kind.

As suggested by Josh in http://github.com/GreenGroup/RMG-Java/issues#issue/10/comment/152512

rwest commented 14 years ago

Not convinced this chebyshev code is safe when it comes to units. What if the user had put the T range in Celsius?

gmagoon commented 14 years ago

Isn't the fundamental problem the low and high points in the ChebyshevPolynomial class? (cf. http://github.com/GreenGroup/RMG-Java/issues/issue/10#issue/10/comment/152485) I think the low and high points should be 800 and 2000 K in Richard's example, but the grid should not be changed.

gmagoon commented 14 years ago

Just to follow up, I think that http://github.com/GreenGroup/RMG-Java/blob/0260e398517e53b766dce9e80738bd53833d8ce1/source/RMG/jing/rxn/PDepRateConstant.java#L271 (the block where TMIN is null) is incorrect, and this is incorrectly giving a Tlow for the Chebyshev of 806.5 K .

rwest commented 14 years ago

I've re-opened this bug as it seems to be unresolved. (The initial problem has been avoided, but by introducing extra problems.)

jwallen commented 14 years ago

A better resolution has been added in af33cca0f19bbe3e590c3e5382723777ce151272 and 48064ce61389a693137e9e5d2e9fa11f2c341ba3. In short, RMG can now automatically distribute the grid of temperature and pressure points for any interpolation model. For Chebyshev it uses the Gauss-Chebyshev grid; for others it distributes the points linearly on 1/T and log P domains. You can also specify the grid manually.

Finally, we now use the interpolation model within RMG, rather than always using bilinear interpolation, to be more consistent and to avoid the temperature out of range error.