Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Out-of-bounds temperature treatment for NASA polynomials #255

Open
mgoodson-cvd opened this issue Feb 29, 2024 · 15 comments · May be fixed by #257
Open

Out-of-bounds temperature treatment for NASA polynomials #255

mgoodson-cvd opened this issue Feb 29, 2024 · 15 comments · May be fixed by #257

Comments

@mgoodson-cvd
Copy link
Collaborator

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:

image

@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
Copy link

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
Copy link
Collaborator

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
Copy link
Collaborator Author

@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
Copy link

Thank you for explaining.

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

Thanks

DL

@mgoodson-cvd
Copy link
Collaborator Author

@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
Copy link

@mgoodson-cvd Thank you for letting me know.

@grgbellasvki
Copy link
Collaborator

@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 added a commit to mgoodson-cvd/Mutationpp that referenced this issue Mar 7, 2024
When the temperature is out-of-bounds with the NASA polynomials,
hold the Cp constant and adjust the enthalpy and entropy accordingly.

Note: this code was mostly authored by Alireza Mazaheri of NASA LaRC.

See mutationpp#255. Closes mutationpp#255.
@mgoodson-cvd
Copy link
Collaborator Author

@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
Copy link

@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 added a commit to mgoodson-cvd/Mutationpp that referenced this issue Mar 7, 2024
Add Matthew Goodson and Alireza Mazaheri to the contributors list.

See mutationpp#255.
mgoodson-cvd added a commit to mgoodson-cvd/Mutationpp that referenced this issue Mar 7, 2024
Update the version to 1.2.0.

See mutationpp#255.
@mgoodson-cvd
Copy link
Collaborator Author

@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
Copy link

@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
Copy link

@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
Copy link
Member

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
Copy link
Collaborator Author

@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
Copy link
Member

@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.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging a pull request may close this issue.

5 participants