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

Implementation of Peng-Robinson EoS #641

Closed
wants to merge 110 commits into from

Conversation

gkogekar
Copy link
Member

@gkogekar gkogekar commented Jun 7, 2019

Changes proposed in this pull request:

  • Implemented Peng-Robinson equation of state
  • Created a new ThermoPhase class - PengRobinsonMFTP
  • Added a new test PengRobinsonMFTP_Test in test/thermo
  • Updated ctml_writer.py to read additional parameter 'acentric_factor' from CTI file

@codecov
Copy link

codecov bot commented Jun 7, 2019

Codecov Report

Merging #641 (64e4acf) into main (81cffde) will increase coverage by 0.10%.
The diff coverage is 71.98%.

❗ Current head 64e4acf differs from pull request most recent head 5a7a938. Consider uploading reports for the commit 5a7a938 to get more accurate results
Impacted file tree graph

@@            Coverage Diff             @@
##             main     #641      +/-   ##
==========================================
+ Coverage   72.39%   72.50%   +0.10%     
==========================================
  Files         364      368       +4     
  Lines       46433    47047     +614     
==========================================
+ Hits        33616    34111     +495     
- Misses      12817    12936     +119     
Impacted Files Coverage Δ
include/cantera/thermo/MixtureFugacityTP.h 0.00% <ø> (ø)
include/cantera/thermo/RedlichKwongMFTP.h 100.00% <ø> (ø)
src/thermo/PengRobinson.cpp 62.54% <62.54%> (ø)
src/thermo/MixtureFugacityTP.cpp 41.75% <67.41%> (+12.98%) ⬆️
src/thermo/RedlichKwongMFTP.cpp 85.92% <84.00%> (+6.73%) ⬆️
include/cantera/thermo/PengRobinson.h 100.00% <100.00%> (ø)
src/thermo/ThermoFactory.cpp 82.45% <100.00%> (+0.05%) ⬆️
test/thermo/PengRobinson_Test.cpp 100.00% <100.00%> (ø)
test/thermo/cubicSolver_Test.cpp 100.00% <100.00%> (ø)
test/thermo/thermoFromYaml.cpp 100.00% <100.00%> (ø)
... and 5 more

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 81cffde...5a7a938. Read the comment docs.

@speth
Copy link
Member

speth commented Jun 8, 2019

Can you please rebase onto the current master and force push, so that this PR does not include a version of the commit "added factories for FlowDevice and Wall objects" which is already merged into the master branch?

Copy link
Member

@decaluwe decaluwe left a comment

Choose a reason for hiding this comment

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

Hi, @gkogekar thanks for adding that additional test. This is generally what I was thinking of. I have a few other suggestions, but it is very very close...

test/thermo/PengRobinsonMFTP_Test.cpp Outdated Show resolved Hide resolved
test/thermo/PengRobinsonMFTP_Test.cpp Outdated Show resolved Hide resolved
test/thermo/PengRobinsonMFTP_Test.cpp Outdated Show resolved Hide resolved
test/thermo/PengRobinsonMFTP_Test.cpp Outdated Show resolved Hide resolved
test/thermo/PengRobinsonMFTP_Test.cpp Outdated Show resolved Hide resolved
Copy link
Member

@bryanwweber bryanwweber left a comment

Choose a reason for hiding this comment

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

Hi @gkogekar and @decaluwe! Thank you so much for this awesome addition, it was clearly a ton of work to put together. Don't mind the number of comments, with this much code to check through, there are always going to be small things. This is overall really good work and I'm looking forward to seeing it in Cantera!

Most of the comments here are about the code style and formatting. I didn't check the implementation/equations thoroughly, but there are a few implementation questions as well. To summarize some of the main comments:

  1. Please use spaces instead of tabs and consistently indent levels by 4 spaces.
  2. Please try to use more descriptive names for variables in some places, the names are kind of ambiguous. Remember that someone else (including you in a few months!) will have to come back and figure out what this code is doing.
  3. I left a comment that a longer documentation should be added to the header file. I'm not sure that's the right place for it though, I wonder if it should go on the website instead?

Otherwise, a few other changes I'd like to see that didn't fit in a line comment:

  • Can this new class be created from Python? I don't see why not, but adding a test to the Python suite would be good.
  • Can this new class be created from a YAML file? Again, don't see why not, but a test of that would be awesome.
  • I'm a little concerned about the code coverage, you're only covering about half of the code you're adding, but we'd like to see closer to 70-80%. I'm especially concerned about the complexity of the cubic solver. Is there any way to add some tests of that solver, even if it is only testing the solver without any physical basis? In other words, is there a way to feed that solver an arbitrary set of coefficients of a cubic equation to test that it works properly?

Please let me know any questions/concerns that you have!

include/cantera/thermo/PengRobinsonMFTP.h Outdated Show resolved Hide resolved
include/cantera/thermo/PengRobinsonMFTP.h Outdated Show resolved Hide resolved
include/cantera/thermo/PengRobinsonMFTP.h Outdated Show resolved Hide resolved
include/cantera/thermo/PengRobinsonMFTP.h Outdated Show resolved Hide resolved
include/cantera/thermo/PengRobinsonMFTP.h Outdated Show resolved Hide resolved
test/data/co2_PR_example.cti Outdated Show resolved Hide resolved
test/data/co2_PR_example.cti Outdated Show resolved Hide resolved
test/thermo/PengRobinsonMFTP_Test.cpp Outdated Show resolved Hide resolved
test/thermo/PengRobinsonMFTP_Test.cpp Outdated Show resolved Hide resolved
test/thermo/PengRobinsonMFTP_Test.cpp Outdated Show resolved Hide resolved
@decaluwe
Copy link
Member

Hi @bryanwweber, thanks for these comments! @gkogekar will have a more detailed reply and fix some of these easier issues. Most are very clear, but I suppose we need to decide, re: documentation, if we want it in header files or on the website.

I think long term the website is the better location, save for a very basic description in the header file. Then the header file can just link to the relevant section on the website. In the short term, this would make the website documentation oddly imbalanced (i.e. a whole lot about Peng-Robinson, but not much else), but I suppose that is unavoidable, and not that big a problem.

So should we just agree that she should add this info to the website?

@bryanwweber
Copy link
Member

@decaluwe Yes, I think the website would be a better place. We do intend to transfer all of the thermo model documentation there eventually, maybe the page that we add could note that most of the documentation is still in C++ and will be moved, to explain why P-R is the only one there 😄

test/thermo/PengRobinsonMFTP_Test.cpp Outdated Show resolved Hide resolved
include/cantera/thermo/PengRobinsonMFTP.h Outdated Show resolved Hide resolved
src/thermo/PengRobinsonMFTP.cpp Outdated Show resolved Hide resolved
include/cantera/thermo/PengRobinsonMFTP.h Outdated Show resolved Hide resolved
src/thermo/PengRobinsonMFTP.cpp Outdated Show resolved Hide resolved
src/thermo/PengRobinsonMFTP.cpp Outdated Show resolved Hide resolved
Copy link
Member

@speth speth left a comment

Choose a reason for hiding this comment

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

This looks like some good work toward implementing the Peng-Robinson equation of state. I have two high-level requests, which I think will help improve this PR.

First, many of the methods appear to be exactly the same as in RedlichKwongMFTP (e.g. enthalpy_mole, entropy_mole, calcDensity, setToEquilState, setTemperature, etc.), and there are others that could be refactored to extract common portions (e.g. the cubic equation solver, the getCoeff method for reading from critProperties.xml). I think it would be worth refactoring to eliminate this redundancy. The only question is whether the common methods should just be moved up to MixtureFugacityTP, or whether there is a need for a new intermediate class. My preference would be for the former.

Second, do we need to have XML-based construction? Since the YAML format is now implemented, can we make this the first class to be part of that format only?

include/cantera/thermo/PengRobinsonMFTP.h Outdated Show resolved Hide resolved
test/thermo/PengRobinsonMFTP_Test.cpp Outdated Show resolved Hide resolved
test/thermo/PengRobinsonMFTP_Test.cpp Outdated Show resolved Hide resolved
@ischoegl
Copy link
Member

ischoegl commented Oct 6, 2019

I agree with @bryanwweber’s comments: thanks for putting in a substantial amount of work!

Beyond, It would be great if this PR could include an introductory example that illustrates usage (python would be ideal). In the spirit of #652, cantera users would benefit from a better overview of what different thermo objects are available and how to use them. Perhaps the example could, as an example, juxtapose ideal gas (default), Redlich-Kwong and the new Peng-Robinson EoS?

@decaluwe decaluwe force-pushed the Peng-Robinson branch 3 times, most recently from 61f04cc to d879455 Compare December 20, 2019 17:08
@decaluwe
Copy link
Member

decaluwe commented Jan 6, 2020

@speth and @bryanwweber - a question about how non-ideal EoS are handled in YAML. I was looking to add an implementation of PengRobinson to the nDodecane_Reitz.yaml file that already ships with cantera. I've noticed that the pureFluid parameters have now moved from the phase definition to the species declaration. I get the logic of this move (and I'm sure @speth asked me about it when making the change), but have a couple of questions:

  1. These are the pure fluid parameters, which go on the diagonal of the multi-component interaction parameter matrix. In the absence of any species interaction parameters, the off-diagonal terms for any pair of species j,k are calculated as a geometric average of the pure fluid parameters (sqrt(a_i*a_j)). However, if a user wants to provide a cross-fluid interaction parameter (for example, if they have solubility data), are they able to do so? Where is the field?
  2. If I wanted to also have a PengRobinson phase declared in this file, can I add a 2nd equation-of-state:model field? Moreover, if I load the ideal-gas phase from this file, is the equation-of-state field simply ignored for the species declarations?

Thanks,
Steven

@bryanwweber
Copy link
Member

bryanwweber commented Jan 7, 2020

For Redlich-Kwong, the cross-fluid parameters are specified in a field called binary-a: https://cantera.org/documentation/dev/sphinx/html/yaml/species.html#redlich-kwong
You can see an example of that here:

- name: H2O
composition: {H: 2, O: 1}
thermo:
model: NASA7
temperature-ranges: [200.0, 1000.0, 3500.0]
data:
- [4.19864056, -0.0020364341, 6.52040211e-06, -5.48797062e-09, 1.77197817e-12,
-30293.7267, -0.849032208]
- [3.03399249, 0.00217691804, -1.64072518e-07, -9.7041987e-11, 1.68200992e-14,
-30004.2971, 4.9667701]
equation-of-state:
model: Redlich-Kwong
a: [1.7458e+08 bar*cm^6/mol^2*K^0.5, -8.0e+04 bar*cm^6/mol^2*K^-0.5]
b: 18.18 cm^3/mol
binary-a:
H2: [4 bar*cm^6/mol^2*K^0.5, 40 bar*cm^6/mol^2*K^-0.5]
CO2: 7.897e7 bar*cm^6/mol^2*K^0.5

That file (test/data/thermo-models.yaml) is a really good reference for some of the more esoteric formats.

Because of the YAML format, you cannot duplicate keys in a mapping, so you can't do something like this:

equation-of-state:
  model: Redlich-Kwong
  model: Peng-Robinson

or

species:
- name: H2O
  equation-of-state:
    model: Redlich-Kwong
  equation-of-state:
    model: Peng-Robinson

So I think you'd need a whole new species definition to add the PR parameters. It seems like the best bet is to make a new file for the PR case. I wonder if it would be better to have a simpler (although obviously less useful) example input file for these cases.

@speth
Copy link
Member

speth commented Jan 7, 2020

2. I wanted to also have a PengRobinson phase declared in this file, can I add a 2nd equation-of-state:model field? Moreover, if I load the ideal-gas phase from this file, is the equation-of-state field simply ignored for the species declarations?

No, you can't have multiple equation-of-state fields for a single species (or duplicate keys in a YAML map in general). Currently, you need to use a different species definition if you want to store different equation of state data.

I hadn't really run across any examples where it seemed useful enough to implement a way of doing this while implementing the YAML format, but I did have one idea in mind that we could consider, which would be to allow the equation-of-state field to be a list, where each item in the list is a map like the current equation-of-state entry, e.g.

equation-of-state: 
- model: Redlich-Kwong 
  a: 1.7458e+08 bar*cm^6/mol^2*K^0.5
  b: 18.18 cm^3/mol 
- model: Peng-Robinson
  a: 2.34e+09
  b: 12.5 cm^3/mol 

Then, when adding the species to the phase, the thermo model could decide which entry from this list to use. I'm not sure how much this ends up complicating the process of adding the species. In this particular case, I think the implementation would be pretty simple. There might be others where it ends up being a mess, but I guess we could always disallow the behavior in those cases.

@decaluwe
Copy link
Member

decaluwe commented Jan 7, 2020

Thanks for the feedback, @bryanwweber and @speth.

My thoughts are not fully formed here, but:

  1. The first thing that occurs to me is that storing EoS-dependent info with the species seems like a significant departure, philosophically, in terms of Cantera's approach to modularity. In this approach, one of these species is no longer a building block that can be used interchangeably to construct any number of phases.

  2. There is certainly a question about how common this use case is (i.e. "Just how frequently will someone want to construct a yaml file with multiple MixtureFugacityTP phases drawing upon the same species?").

I would guess that the answer to the latter is "not very often," but I also think it is a nice goal to maintain the original spirit of the software, if the implementation of equation-of-state as a list is not too big a mess. In the cases where someone does want to implement multiple non ideal EoS, requiring them to duplicate all species seems like a pain.

@bryanwweber - I could certainly construct a second yaml file for just this phase, but I'm thinking of making an example that can demonstrate, more generally, Cantera's non ideal EoS capabilities. Having a single input file for this purpose would be nice, both so that users have a one-stop reference, and to avoid "input file bloat." if possible. 😁

Thanks again, to you both!

@speth
Copy link
Member

speth commented Jan 7, 2020

  1. The first thing that occurs to me is that storing EoS-dependent info with the species seems like a significant departure, philosophically, in terms of Cantera's approach to modularity. In this approach, one of these species is no longer a building block that can be used interchangeably to construct any number of phases.

While it is a departure from the XML/CTI format, and setting aside the current limitations of the equation-of-state data structure, I think storing such information with the species improves reusability / modularity. A few examples:

  • I want to add a species to an R-K phase. In the CTI/XML approach, I have to provide both the species definition and add new parameters to the phase definition as well.
  • I want to create multiple R-K phases that includes different subsets of species. In the CTI/XML approach, I would have to duplicate all of the R-K parameters in each phase definition.
  • I want to create phases using different equations of state (i.e. R-K and P-R) based on more fundamental species data (i.e. critical temperature and pressure). In the CTI/XML approach, I would have to duplicate this data in each phase definition.

I think there is also a parallel with how we handle species transport data - the parameters used by multiple transport models are part of the species definitions. The same limitation and possible extension of allowing this to be a list with data for multiple models applies here as well, for example if we were to introduce a transport model that didn't want the Lennard-Jones parameters.

@decaluwe
Copy link
Member

decaluwe commented Jan 7, 2020

These are good points. Regardless of whether we enable listing of multiple EoS for a given species, I will proceed at the moment with a separate YAML for the Peng–Robinson tests and examples.

Would the preference be that I eliminate all cti files from the PR? I'm guessing yes.

@decaluwe
Copy link
Member

decaluwe commented Jan 8, 2020

note the typo in my last post has been corrected. I've added a YAML implementation. Should I remove all cti and xml files and references?

@speth
Copy link
Member

speth commented Jan 9, 2020

note the typo in my last post has been corrected. I've added a YAML implementation. Should I remove all cti and xml files and references?

Yes, I think any new models introduced at this point should use only the YAML format.

@decaluwe
Copy link
Member

Hi all, just wanted to get everyone on the same page, for this PR.

By my reckoning, there are four main changes required to get this PR on solid footing for merging. The first two relate more generally to a wish for greater test coverage:

  1. A test of the general performance of the cubic solver. @gkogekar is working on this, at present, to be added to the C++ test suite.
  2. Some tests added to the python test suite. I am working on this.
  3. @ischoegl requested an example of the class performance and capabilities. My current thinking is to add this to the Jupyter notebook example, showing more generally the non-ideal capabilities (both the pending PengRobinson and RedlichKwongMFTP). This has the side effect, of course, of making it a differnt PR.
  4. There are still a few common functions in PengRobinson and RedlichKwongMFTP which could be abstracted up to the MixtureFugacityTP class. These remaining functions all involve setTemperature, which for some reason fails the PengRobinson test suite when I move it. I have this commit pushed to my repo: decaluwe@bbedc50

Thoughts on these?

@speth and @bryanwweber it would also be good to know the time frame on the next release. Obviously I'd like to get this merged before then.

Thanks!

@ischoegl
Copy link
Member

ischoegl commented Jan 16, 2020

@decaluwe

@ischoegl requested an example of the class performance and capabilities. My current thinking is to add this to the Jupyter notebook example, showing more generally the non-ideal capabilities (both the pending PengRobinson and RedlichKwongMFTP). This has the side effect, of course, of making it a differnt PR.

While I like the idea of using jupyter in general, I believe it may be illustrative to have a really simple/bare-bones illustration for how EoS behavior changes for the different thermo models. From what I've seen, there's currently only a NonIdealShockTube.ipynb, which throws in several more complications. I was hoping to see something more along an intro to thermo perspective instead of a combustion application. There's currently no such thing in the thermo folder (or the Python examples of the Cantera package).

PS: Looks like my comment isn't too far from what's already suggested. My main point is to make it a new example in a different folder.

throw CanteraError("PengRobinson::setSpeciesCoeffs",
"Unknown species '{}'.", species);
}

Copy link

Choose a reason for hiding this comment

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

According to the definition proposed by Pitzer, acentric factor is not likely to be less than -1. So is it more appropriate to detect the effectiveness of ω and report the incorrect input here?

Copy link
Member

Choose a reason for hiding this comment

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

Hi, @Yuji-1 - thanks for this comment!

There is always a balance in what the user is allowed to input, in terms of what is known/thought to be "correct," and what a user might want to try out, for some reason perhaps not known to us. Think of it as a balance between "rules" and "freedom."

In this case, so long as an acentric factor less than -1 does not lead to errors or failures in the code, my preference would be to allow it.

Copy link
Member

@bryanwweber bryanwweber left a comment

Choose a reason for hiding this comment

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

A few unaddressed comments here, plus one minor formatting comment. Thanks @gkogekar!

src/thermo/PengRobinson.cpp Outdated Show resolved Hide resolved
Comment on lines 488 to 493
double a0 = 0, a1 = 0;
if (item2.second.isScalar()) {
a0 = units.convert(item2.second, "Pa*m^6/kmol^2");
}
setBinaryCoeffs(item.first, item2.first, a0, a1);
}
Copy link
Member

Choose a reason for hiding this comment

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

As @speth said here, you're declaring a1 to be zero and then passing it into setBinaryCoeffs right away, but is it supposed to be something else?

src/thermo/PengRobinson.cpp Outdated Show resolved Hide resolved
@speth
Copy link
Member

speth commented Apr 20, 2021

I put a good bit of effort into cleaning up the commit history on this branch to eliminate some of the unnecessary churn (commits switching line endings from LF to CRLF and back, commits that mistakenly added lines containing merge conflict markers, etc.) and pushed that branch to my fork (speth/Peng-Robinson). I also added several new tests to cover some parts of the code that were not being exercised, and modified the conditions of some of the other tests to look at situations involving multiple species or where the non-ideal behavior should be more significant. As a result, there are now 6 test cases that do not pass:

  • PengRobinson_Test.chemPotentials_RT
  • PengRobinson_Test.activityCoeffs
  • PengRobinson_Test.cpValidate
  • PengRobinson_Test.partialMolarCpFiniteDifference
  • PengRobinson_Test.partialMolarPropertyIdentities
  • PengRobinson.lookupSpeciesPropertiesMissing

@gkogekar, my suggestion at this point would be to reset your branch to match mine, and close this PR. Once you've had a chance to resolve the issues leading to these test failures, you can then open a new PR and we won't have this enormous comment history bogging things down.

@gkogekar
Copy link
Member Author

Closing this PR and creating a new one.

@gkogekar
Copy link
Member Author

A revised PR is created here: #1047

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.

7 participants