Out-of-bounds temperature treatment for NASA polynomials

[mgoodson-cvd]

The current behavior with the NASA thermodynamic polynomials is to use them directly, even if the temperature is out-of-bounds. This can be extremely problematic as the polynomials go wild outside of the bounds. (This issue was already partially raised in #117). Note that each NASA-9 table entry provides its own temperature bounds, though most of the McBride data is 200 K - 6000 K (or 20,000 K).

The common treatment for out-of-bounds temperature with NASA polynomials is to hold Cp constant at the last in-bounds temperature, e.g., from the FUN3D version 13.4 manual:

@AlirezaMazaheri provided me with his initial implementation of this in Mutation++ and his permission to share it. I have updated it to work with the latest version of Mutation++ and am currently testing it. I hope to be able to push it soon.

[domilanza2002]

Good evening, I am trying to shift the reference temperature from 25°C to 0°C. In particular, I am using python, so the branch pypi. The formula I am using is: h_new(T,P)=h(T,P)+(h-h0)(T,P) This works for RRHO and NASA-7, but when I try to use it with NASA-9, I get a huge (h-h0) enthalpy. In particular I am doing: import mutationpp as mpp mix=mpp.Mixture("nitrogen2") mix.equilibrate(298.15,101325) print(mix.mixtureHMinusH0Mass()) I get the value: 65604568531370.125

Do you think this is related to the issue you are describing? In that case, could you share the part of the code that must be edited?

Thanks

[grgbellasvki]

Hi @mgoodson-cvd. That would be an interesting addition. I can review the merge request. I wonder whether adding an option to choose the default behavior is important, to avoid breaking other codes. Even though yours is the most natural choice, it is not what people always do. At VKI, we had to add fixes like the one you mention, but they were always in the client code.

@domilanza2002 Probably these two issues are indeed related.

[mgoodson-cvd]

@domilanza2002 That is precisely the issue, because mixtureHMinusH0Mass gets the mixture enthalpy at approximately 0 K, which for NASA-9 polynomials is going to be wildly incorrect without the proposed fix.

@grgbellasvki The switch is an interesting idea, though I can't imagine ever running without it on. I would certainly want this to be the default behavior to avoid issues like @domilanza2002. What other choices for out-of-bounds treatment are there?

[domilanza2002]

Thank you for explaining.

When do you think you will be able to push the new version?

Thanks

DL

[mgoodson-cvd]

@domilanza2002 I don't have a timetable at this point. It will depend on how @grgbellasvki and others want to proceed. I will try to push at least a working version to my fork in the next week or so, though with no tests or options to turn off.

[domilanza2002]

@mgoodson-cvd Thank you for letting me know.

[grgbellasvki]

@mgoodson-cvd I was thinking more of throwing an error if the temperature is out-of-bounds as the default behavior and the switch would enable extrapolation. I agree that what you recommend is a significant improvement compared to what it is now in the code.

If you push a merge request, I can review the code and try a few LTE and B' calculations which are of interest for many users.

As far as alternatives, it really depends on the application. I know for example a team which freezes enthalpy below a certain value to assist chemistry thermodynamics. Another approach would be extrapolating at low Ts using RR. Warning the user to provide more data is often the best and safest choice.

Thanks!

[mgoodson-cvd]

@grgbellasvki @domilanza2002 I've opened #257 with an initial draft of the code.

@grgbellasvki I don't think erroring out when out-of-bounds would be desired. For one thing, some elements like electrons have a lower temperature bound of 298.15 K; yet the Cp is already constant for electrons so there really isn't a "lower" bound. Also, the polynomials aren't that bad for some range of temperatures outside of the bounds; e.g., for air, NASA-9 is still sane at 160 K, but then gets bad very quickly as the polynomials go haywire.

In regards to warning the user, it sounds like a good idea but in practice could get extremely verbose, such as for the electron case mentioned above.

I guess the best path forward would be to add a flag to the mixture inputs to control the out-of-bounds behavior, e.g.,

<mixture thermo_db="NASA-9" nasa_oob="constant_cp">

where the two current options for nasa_oob are constant_cp and none, where none is the old behavior. I would want to switch the default behavior to constant_cp though, so we would want to consider making this at least a minor version increment.

[domilanza2002]

@mgoodson-cvd Thank you! Since I am using the python interface, to use the new version what should I do? Should I change the file in the src folder and build again?

Thank you and sorry for the bothering

[mgoodson-cvd]

@domilanza2002 Yes, you will have to rebuild using my branch, https://github.com/mgoodson-cvd/Mutationpp/tree/255-out-of-bounds-treatment-for-nasa-polynomials

You have a few options. You can just manually update the NasaDB.h file in your repo and replace it with mine. Or you can add my fork as a remote and checkout my branch. Or you could download my fork directly as a new folder.

[domilanza2002]

@mgoodson-cvd Thank you! If you're doing a commit about this I might suggest you add the function to shift the reference temperature too(issue 69). It should be just a few lines of code.

Thanks

[AlirezaMazaheri]

@mgoodson-cvd If you are adding an option I would suggest to use one that is more descriptive than constant_cp, because it could be misinterpreted. How about extrapolate_out_of_range_T? Along the comment made by @grgbellasvki, I think the none option should return an error when the database is accessed beyond the valid T range.

[jbscoggi]

Hey All, just to add my two cents. I agree that cp clipping is probably the "right" solution and I also agree that having an option for the user is also desired. I see three possibilities that would be useful:

1. clip_cp_at_temperature_bounds (default)
2. extrapolate
3. raise_error

Option 1 is the cp clipping described by @mgoodson-cvd. Option 2 is the current behavior. Option 3 is to explicitly crash by raising an exception. I would call this option "when_outside_temperature_range" in the mixture file. Putting this together with any of the three options is pretty clear IMO. Ideally, this could be implemented with a strategy pattern that doesn't require a bunch of if statements every time you hit this decision point.

[mgoodson-cvd]

@AlirezaMazaheri @jbscoggi Thanks both of y'all for the feedback here. I will give this a shot when I get a chance. I was already thinking about how to handle the "if" statements, either through templating or multiple versions of the NasaDB class. I'll let you know when I have something working.

[jbscoggi]

@mgoodson-cvd, I think there are probably two approaches that work well. One is using the decorator pattern in which you have the NasaDB class as it is (this represents the extrapolate option from above. The other options are then implemented as decorators, which are classes that inherit from NasaDB and own a NasaDB type and wrap all the relevant function calls by first checking the temperature range and then if deciding how to call the vanilla NasaDB type to get the result. Depending on the user-selected option, you would create the appropriate object at runtime when loading the database. Another option is the strategy pattern in which the NasaDB type would be modified to take an ExtrapolateStrategy type (or whatever you want to call it) that implements some interface that allows you to abstract out the details of how the database is extrapolated. Not clear what that interface should look like off the top of my head, but that's the general idea. Again, at start-up you would create the appropriate strategy object and pass into the database type. The best approach might actually be some combination of those ideas, where you have a delegate class that is responsible for checking the temperature range and using a given strategy to implement the details of what it actually does when the temperature is outside of the allowable range. This has the nice feature of being easily extensible to new strategies without having to modify the existing database class.