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

Individual Reactor Initialization - Segmentation Fault Fix. #1063

Merged
merged 8 commits into from
Jun 26, 2021

Conversation

anthony-walker
Copy link
Contributor

There was a segmentation fault when calling initialize on an individual reactor, this is a fix for that by setting it to the specified time if a ReactorNet has not been provided or used.

I specifically changed ReactorNet::updateConnected in both the .cpp and .h files to accept a time arguement. If the reactor network is not available, it uses the specified or default time value of the reactor specific initialize function.

Checklist

  • The pull request includes a clear description of this code change
  • Commit messages have short titles and reference relevant issues
  • Build passes (scons build & scons test) and unit tests address code coverage
  • Style & formatting of contributed code follows contributing guidelines
  • The pull request is ready for review

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.

Thanks for the quick fix @anthony-walker! Is it necessary to pass in t0 to updateConnected? If a reactor is not connected to a Network, any times greater than zero aren't really meaningful, and if it is connected to a network this won't be a problem. I suggested just setting the value to zero in the change below.

Can you please also add a test for this, so we don't re-introduce the bug by accident, and we can make sure that it's fixed?

@@ -187,7 +187,7 @@ void Reactor::updateConnected(bool updatePressure) {
m_thermo->saveState(m_state);

// Update the mass flow rate of connected flow devices
double time = m_net->time();
double time = (m_net!=NULL) ? m_net->time() : t0;
Copy link
Member

Choose a reason for hiding this comment

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

Suggested change
double time = (m_net!=NULL) ? m_net->time() : t0;
double time = (m_net != NULL) ? m_net->time() : 0.0;

Cantera's style is to have space around the operators.

Copy link
Member

Choose a reason for hiding this comment

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

Also, you should compare to nullptr rather than NULL to be more in line with C++11.

@bryanwweber
Copy link
Member

Actually, thinking about it a little more... The other option here is to throw an error if the network is not connected. Did you consider that option? If you did, and you rejected it, can you explain why? Thanks!

@anthony-walker
Copy link
Contributor Author

@bryanwweber I considered both using strictly zero and throwing an error. I chose to pass the argument along because it is presently used in all of the initialize functions. I suspect you are correct that it will be zero most of the time as if folks are integrating they will use a network. I rejected throwing the error because it will be useful in my case for writing some unit tests where I strictly need values from the reactor. I also believe that in lieu of throwing an error it should just be a protected function and ReactorNet can be made a friend of it to allow it to be called in the conventional way.

I will also add a test for which ever direction you feel we should pursue.

@codecov
Copy link

codecov bot commented Jun 17, 2021

Codecov Report

Merging #1063 (9fec84f) into main (9de22e8) will increase coverage by 0.14%.
The diff coverage is 100.00%.

Impacted file tree graph

@@            Coverage Diff             @@
##             main    #1063      +/-   ##
==========================================
+ Coverage   72.92%   73.07%   +0.14%     
==========================================
  Files         357      358       +1     
  Lines       47120    46956     -164     
==========================================
- Hits        34363    34311      -52     
+ Misses      12757    12645     -112     
Impacted Files Coverage Δ
src/zeroD/Reactor.cpp 84.81% <100.00%> (+1.11%) ⬆️
test/zeroD/test_zeroD.cpp 100.00% <100.00%> (ø)
src/base/global.cpp 72.38% <0.00%> (-2.86%) ⬇️
src/thermo/StoichSubstance.cpp 83.69% <0.00%> (-2.18%) ⬇️
src/base/units.h 92.10% <0.00%> (-0.88%) ⬇️
test_problems/cathermo/stoichSub/stoichSub.cpp 89.13% <0.00%> (-0.46%) ⬇️
src/thermo/SpeciesThermoFactory.cpp 78.81% <0.00%> (-0.38%) ⬇️
test_problems/diamondSurf/runDiamond.cpp 94.80% <0.00%> (-0.32%) ⬇️
src/base/xml.cpp 72.20% <0.00%> (-0.19%) ⬇️
src/thermo/HMWSoln.cpp 81.07% <0.00%> (-0.18%) ⬇️
... and 24 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 9de22e8...9fec84f. Read the comment docs.

Copy link
Member

@ischoegl ischoegl left a comment

Choose a reason for hiding this comment

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

I think the proposed fix looks good other than the usual formatting quibbles.

Edit: 👍 on @bryanwweber ’s request to add a unit test. Perhaps you can add the one you already started? (Minus stuff that is specific to the adaptive preconditioner).

@@ -88,7 +88,7 @@ void Reactor::initialize(doublereal t0)
m_thermo->restoreState(m_state);
m_sdot.resize(m_nsp, 0.0);
m_wdot.resize(m_nsp, 0.0);
updateConnected(true);
updateConnected(true,t0);
Copy link
Member

Choose a reason for hiding this comment

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

Another missing space.

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.

Seems like a pretty straightforward fix.

One suggestion for the test case would be to make sure that the reactor still works correctly if it's initialized without a ReactorNet and then later added to one. I don't see a reason why this wouldn't work, but it would be good to check.

@@ -193,7 +193,8 @@ class Reactor : public ReactorBase
//! `true` for reactors where the pressure is a dependent property,
//! calculated from the state, and `false` when the pressure is constant
//! or an independent variable.
virtual void updateConnected(bool updatePressure);
//! @param t0 initialization time for the reactor if not obtained from the network
virtual void updateConnected(bool updatePressure, double t0 = 0.0);
Copy link
Member

Choose a reason for hiding this comment

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

No spaces around the = when it's defining a default value.

@@ -187,7 +187,7 @@ void Reactor::updateConnected(bool updatePressure) {
m_thermo->saveState(m_state);

// Update the mass flow rate of connected flow devices
double time = m_net->time();
double time = (m_net!=NULL) ? m_net->time() : t0;
Copy link
Member

Choose a reason for hiding this comment

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

Also, you should compare to nullptr rather than NULL to be more in line with C++11.

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.

Thanks @anthony-walker! These changes look good to me now, with the exception of the comment below about the test.

test/zeroD/test_zeroD.cpp Outdated Show resolved Hide resolved
Copy link
Member

@ischoegl ischoegl left a comment

Choose a reason for hiding this comment

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

One more quibble … your commit messages don’t follow guidelines (per CONTRIBUTING.md):

Use descriptive commit messages (summary line of no more than 72 characters, followed by a blank line and a more detailed summary, if any)

Would you mind updating those messages? Sorry about being nitpicky … it’s simple to fix using ‘reword’ when rebasing.

There was a segmentation fault when calling initialize on an individual reactor, this is a fix for that by setting it to the specified time if the network is not available
I adapted the example ./interfaces/cython/cantera/test_reactor.py::test_equilibrium_HP to C++ and performed the reactor initialization with an arbitrary value to show that it no longer creates a segmentation fault and also still integrates successfully in a network post individual initialization. I added this to it's own folder called zeroD because I could not find C++ instances of Reactor tests elsewhere.
This commit changes the former test to use Reactor classes in lieu of IdealGasConstPressureReactor classes because it is the base class. It also changes some of the comparisons between the two reactor objects.
Copy link
Member

@ischoegl ischoegl left a comment

Choose a reason for hiding this comment

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

Quick comment: the new test failed on the windows runners and needs tweaks, see

cl /Fobuild\test\zeroD\test_zeroD.obj /c test\zeroD\test_zeroD.cpp /EHsc /MD /nologo /D_SCL_SECURE_NO_WARNINGS /D_CRT_SECURE_NO_WARNINGS /W3 /DNDEBUG /Iext\googletest\googletest\include /Iext\googletest\googlemock\include /Iinclude /I3rdparty\boost test_zeroD.cpp
test\zeroD\test_zeroD.cpp(44): error C2131: expression did not evaluate to a constant
test\zeroD\test_zeroD.cpp(44): note: failure was caused by call of undefined function or one not declared 'constexpr'
test\zeroD\test_zeroD.cpp(44): note: see usage of 'Cantera::Reactor::neq'
test\zeroD\test_zeroD.cpp(45): error C2131: expression did not evaluate to a constant
test\zeroD\test_zeroD.cpp(45): note: failure was caused by call of undefined function or one not declared 'constexpr'
test\zeroD\test_zeroD.cpp(45): note: see usage of 'Cantera::Reactor::neq'
scons: *** [build\test\zeroD\test_zeroD.obj] Error 2
scons: building terminated because of errors.

Specifically, it ends lines of at 88 characters and adds space after inline comments. It also changed the allocation of arrays to satisfy CI.
Copy link
Member

@ischoegl ischoegl left a comment

Choose a reason for hiding this comment

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

Thanks for the CI fix suggestion. While CI has not finished at the moment, here's a comment.

test/zeroD/test_zeroD.cpp Outdated Show resolved Hide resolved
ischoegl
ischoegl previously approved these changes Jun 20, 2021
double tol = 1e-14;
EXPECT_NEAR(state1[0], state2[0], tol);
EXPECT_NEAR(state1[1], state2[1], tol);
for(size_t i = 3; i < reactor1.neq(); i++)
Copy link
Member

Choose a reason for hiding this comment

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

One thing I hadn't noticed before: why are you skipping i=2? Shouldn't all state variables be comparable?

Also, as minor formatting suggestions: usually Cantera doesn't use a line break before { in a for loop (they're typically on the same line). There definitely should be a space after the for, though. Beyond, it may also help to add a blank line before each of your comments, as it makes your code easier to parse visually.

I'll let @bryanwweber make the final approval, as he was the one who started the review.

Copy link
Contributor Author

@anthony-walker anthony-walker Jun 25, 2021

Choose a reason for hiding this comment

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

So, I adapted this example from test_equilibrium_HP where pressure is held as a constant and not compared. When it is compared it causes the test to fail based on the specified tolerance in the test. I am not entirely certain why but I assumed is was due to the methods used and was justified at the time. This example is the same in that U is held constant but varies slightly on the order of 1e-8. I will look further into it though for a more thorough explanation.

Copy link
Member

Choose a reason for hiding this comment

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

It’s actually interesting that this unit test works to the specified tolerances at all: the test compares a result that involves reaction kinetics with one that uses equilibrium. While there shouldn’t be much of a difference for these simple mixtures in theory, it’s not completely intuitive that comparisons are that close. I’m saying this as I’ve seen significant differences between equilibrium results and kinetics calculations for some more complex mixtures.

On the other hand, I f this is the exact equivalent to a Python example, why don’t you just use the ‘hp’ equilibrium case, plus IdealGasConstPressureReactor objects for this unit test? That it works for HP and not for UV could just be a quirk of some solvers.

Copy link
Contributor Author

@anthony-walker anthony-walker Jun 26, 2021

Choose a reason for hiding this comment

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

@ischoegl I could reduce the tolerance and the unit test would work. I believe it passes for 1e-7 but I didn't think that was a good answer. It doesn't work for HP either though if that was unclear. Also, Reactor was selected to be the test reactor because it is the base class and not a derived class.

Copy link
Member

Choose a reason for hiding this comment

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

IMO, 1e-7 is a pretty good tolerance on the in internal energy... We're talking about a value on the order of 1e6 right?

Copy link
Member

@ischoegl ischoegl Jun 26, 2021

Choose a reason for hiding this comment

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

It doesn't work for HP either though if that was unclear.

no, I guess it wasn’t clear to me … in that case, I’d agree with @bryanwweber and just reduce the tolerance. 1e-14 is really tight if you’re not comparing apples to apples.

Otherwise, please have another look at the formatting as suggested. Also, note that there are some additional white spaces missing on line 17, I.e. double P0 = 10*OneAtm; should be double P0 = 10 * OneAtm;

Copy link
Contributor Author

Choose a reason for hiding this comment

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

@ischoegl I will fix the formatting. I always forget that. Hopefully I will catch on to that eventually. I will also make the change to 1e-7. @bryanwweber the values are pretty large, 971244.30 to be more specific.

Copy link
Member

Choose a reason for hiding this comment

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

@anthony-walker Right, so that's ≈1e6. A tolerance of 1e-14 is on the order of machine precision, and applying several arithmetic operations (integrating or solving equilibrium) is going to decrease the precision of the answer. Since there's different methods being compared here, a tolerance down at machine precision is too tight. On top of that, since the values are so large, at some point the number of digits is no longer physically meaningful. If you look at some of the tests in the regression tests, the tolerances are on the order of 1e-2 to account for differences between compilers, platforms, etc.

Anyhow, all of which is to say that sometimes, changing the tolerance is a legitimate answer 😊

Copy link
Contributor Author

Choose a reason for hiding this comment

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

@bryanwweber Thank you for the explanation, it makes more sense now. I will have to look through those tests. Have all concerns been addressed with the PR then?

Copy link
Member

Choose a reason for hiding this comment

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

It’s actually interesting that this unit test works to the specified tolerances at all: the test compares a result that involves reaction kinetics with one that uses equilibrium. While there shouldn’t be much of a difference for these simple mixtures in theory, it’s not completely intuitive that comparisons are that close. I’m saying this as I’ve seen significant differences between equilibrium results and kinetics calculations for some more complex mixtures.

Since Cantera computes reverse reaction rates from equilibrium constants, if you have a reaction mechanism with at least as many independent, reversible reactions as species, I think you should always get "exactly" the same result from both a kinetic simulation out to steady state and the equilibrium solver. Given that both the equilibrium and ODE solvers are only accurate to within certain bounds, the level of exactness depends on the tolerances of these solvers, not to something approximating accumulated floating point roundoff.

So, I adapted this example from test_equilibrium_HP where pressure is held as a constant and not compared. When it is compared it causes the test to fail based on the specified tolerance in the test.

For future reference, note that the Python assertNear function is applying a combined relative and absolute tolerance. The GTest ASSERT_NEAR is just an absolute tolerance. To get a relative tolerance, which is generally more useful, we often end up writing something like ASSERT_NEAR(x0, x1, 1e-7 * std::abs(x0)).

test/zeroD/test_zeroD.cpp Outdated Show resolved Hide resolved
It specifically changes the application of the phase and kinetics
objects to the reactor by using "insert" instead of setKineticMgr and
setThermoMgr. It also changes the initial pressure to use a built in
constant in cantera.
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.

Looks good to me, thanks @anthony-walker

@bryanwweber bryanwweber merged commit 44f54d2 into Cantera:main Jun 26, 2021
@anthony-walker anthony-walker deleted the bugfix branch June 29, 2021 20:19
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.

4 participants