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

Remove E0 from Davis Products EOS #372

Merged
merged 18 commits into from
May 14, 2024
Merged

Remove E0 from Davis Products EOS #372

merged 18 commits into from
May 14, 2024

Conversation

jhp-lanl
Copy link
Collaborator

@jhp-lanl jhp-lanl commented May 10, 2024

PR Summary

We were taking the E0 parameter in the Davis Products EOS but not actually using it anywhere.

This MR simply subtracts the E0 from the normal reference isentrope, which I believe should be sufficient and consistent in the rest of the EOS since es is involved in all the energy expressions.

I also added a comment basically repeating the need described in #260 for an iterative procedure to find the CJ energy from a reference state.

NOTE: This is really just a quick fix so I haven't added any unit tests. I hope to add a more complete set of tests later with a bit more of a revamp to the EOS.

PR Checklist

  • Adds a test for any bugs fixed. Adds tests for new features.
  • Format your changes by using the make format command after configuring with cmake.
  • Document any new features, update documentation for changes made.
  • Make sure the copyright notice on any files you modified is up to date.
  • After creating a pull request, note it in the CHANGELOG.md file.
  • LANL employees: make sure tests pass both on the github CI and on the Darwin CI

If preparing for a new release, in addition please check the following:

  • Update the version in cmake.
  • Move the changes in the CHANGELOG.md file under a new header for the new release, and reset the categories.

@jhp-lanl jhp-lanl requested review from chadmeyer and jonahm-LANL May 10, 2024 01:48
@jhp-lanl
Copy link
Collaborator Author

Updated the documentation to reflect the change and explain how to calculate E0 from the CJ state

Copy link
Collaborator

@Yurlungur Yurlungur left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Looks like there's still some TODOs in here. And I'd like @chadmeyer 's eyes on it. We should maybe also add a unit test, and trigger it on Darwin.

That said I'm going to approve now, in case you need to get this in quickly for EAP integration work.

@@ -292,6 +292,7 @@ class DavisProducts : public EosBase<DavisProducts> {
inline void Finalize() {}
static std::string EosType() { return std::string("DavisProducts"); }
static std::string EosPyType() { return EosType(); }
// TODO (JHP): Create helper function to find the CJ state given the reference state
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

👍 I'm in favor of the helper being static too.

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The member function will need to know the input parameters, so if I make it static it will need to take the appropriate parameters in addition to the reference pressure and density.

I'll move this debate to #260

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Given that we now use the shift to set the energy of the CJ state, I am doubly convinced we should be providing a method to compute it, so that the shift can be straightforwardly generated. I suppose it doesn't need to be static necessarily... since you could generate an EOS, compute the energy, then shift it.

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yeah... since the modifier takes an initialized object, it seems most prudent to have it not be static so that you don't have to copy/paste parameters and risk an error.

We might want to discuss what this should look like on the Fortran side though. I'll make a comment in the issue.

@jhp-lanl
Copy link
Collaborator Author

Looks like there's still some TODOs in here. And I'd like @chadmeyer 's eyes on it. We should maybe also add a unit test, and trigger it on Darwin.

I intend to address the TODO in a future MR. While it won't be overly complicated to implement the CJ state solver, it is a bit more work than I need to get things working in EAP. Philosophically, it's also different from the changes here; this is a fix since E0 was irrelevant before, whereas the CJ state solver is a new feature.

@chadmeyer this should be pretty fast to review

@chadmeyer
Copy link
Collaborator

There is a potential issue here. I was talking with Andrew Henrick, and for consistency, there might need to be other changes. For instance, $P_s = \frac{e_s}{V}(k-1+F(V))$. The point here is that if we shift the isentrope, it might mess with the pressure in unexpected ways.

@jhp-lanl
Copy link
Collaborator Author

There is a potential issue here. I was talking with Andrew Henrick, and for consistency, there might need to be other changes. For instance, Ps=esV(k−1+F(V)). The point here is that if we shift the isentrope, it might mess with the pressure in unexpected ways.

We don't use that functional form in the code as it exists right now. We calcualte the isentrope pressure as

PORTABLE_INLINE_FUNCTION Real Ps(const Real rho) const {
    const Real vvc = 1 / (rho * _vc);
    return _pc * std::pow(0.5 * (std::pow(vvc, _n) + std::pow(vvc, -_n)), _a / _n) /
           std::pow(vvc, _k + _a) * (_k - 1.0 + F(rho)) / (_k - 1.0 + _a);
  }

I haven't gone through the math to prove to myself that this is just the derivative of the energy along the isentrope, but it should be right? And if I'm just adding a constant offset to the energy, then the derivative shouldn't change.

@chadmeyer
Copy link
Collaborator

I haven't gone through the math to prove to myself that this is just the derivative of the energy along the isentrope, but it should be right? And if I'm just adding a constant offset to the energy, then the derivative shouldn't change.

In principle I think this is all right. I just went through the whole code again, and it looks like this will be right. But that begs the question of whether we should have an E0 at all here... If we are correct, then it's just an EOS shift, and we already have machinery for that and maybe we just want to remove the option of doing it custom for this EOS?

@jhp-lanl
Copy link
Collaborator Author

I haven't gone through the math to prove to myself that this is just the derivative of the energy along the isentrope, but it should be right? And if I'm just adding a constant offset to the energy, then the derivative shouldn't change.

In principle I think this is all right. I just went through the whole code again, and it looks like this will be right. But that begs the question of whether we should have an E0 at all here... If we are correct, then it's just an EOS shift, and we already have machinery for that and maybe we just want to remove the option of doing it custom for this EOS?

Ah interesting point... that might be the better solution. I'll see if I can get that working in my tests

@jhp-lanl
Copy link
Collaborator Author

I got the energy shift working in my tests and I think that's the most clear way to make the shift happen. I've removed the E0 from the Davis Products EOS and put in the documentation that the user should explicitly use our energy shift modifier to accomplish the task of properly connect the products EOS to the desired reference state.

@chadmeyer
Copy link
Collaborator

I've removed the E0 from the Davis Products EOS and put in the documentation that the user should explicitly use our energy shift modifier to accomplish the task of properly connect the products EOS to the desired reference state.

I approve of this improvement to the implementation!

@@ -292,6 +292,7 @@ class DavisProducts : public EosBase<DavisProducts> {
inline void Finalize() {}
static std::string EosType() { return std::string("DavisProducts"); }
static std::string EosPyType() { return EosType(); }
// TODO (JHP): Create helper function to find the CJ state given the reference state
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Given that we now use the shift to set the energy of the CJ state, I am doubly convinced we should be providing a method to compute it, so that the shift can be straightforwardly generated. I suppose it doesn't need to be static necessarily... since you could generate an EOS, compute the energy, then shift it.

@jhp-lanl
Copy link
Collaborator Author

Merging now that I stomped out all of the instances of E0 in the python bindings

@jhp-lanl jhp-lanl changed the title Utilize E0 in Davis Products EOS Remove E0 from Davis Products EOS May 14, 2024
@jhp-lanl jhp-lanl merged commit a213659 into main May 14, 2024
5 checks passed
@jhp-lanl jhp-lanl deleted the jhp/DavisEOS branch May 14, 2024 02:17
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 this pull request may close these issues.

3 participants