Adding adaptive preconditioner functionality to Cantera by anthony-walker · Pull Request #1010 · Cantera/cantera

anthony-walker · 2021-04-09T15:43:44Z

This code is adding preconditioning capabilities to Cantera/CVODEs solver interface for reactor networks. It is focused on a specific type of preconditioning but the code is being written with extensibility in mind. It adds Preconditioning capabilities to the numerics portion of cantera which contains a "PreconditionerBase" that is used to write future preconditioners. It also contains a class derived from this called "AdaptivePreconditioner" which is the aforementioned specific preconditioner.

There are also a couple of functions added to specific reactors. IdealGasReactor is one example. Currently, I added a function to evaluate the energy equation which is also done inside of evalEqs but I wanted to avoid touching core Cantera code as much as possible. In the future, I would like to collaborate with the appropriate people to stream line this.

Checklist

There is a clear use-case for this code change
The commit message has a short title & references relevant issues
Build passes (scons build & scons test) and unit tests address code coverage
The pull request is ready for review

anthony-walker · 2021-04-09T15:45:36Z

@bryanwweber and @ischoegl, I have updated the pull request to be from my "ap-dev" branch. I will use this in the future.

ischoegl · 2021-04-09T15:50:29Z

Tagging #951.

codecov · 2021-04-09T16:44:16Z

Codecov Report

Merging #1010 (394bc7f) into main (d60e59e) will increase coverage by 0.08%.
The diff coverage is 71.04%.

❗ Current head 394bc7f differs from pull request most recent head e74ffd3. Consider uploading reports for the commit e74ffd3 to get more accurate results

@@            Coverage Diff             @@
##             main    #1010      +/-   ##
==========================================
+ Coverage   67.93%   68.01%   +0.08%     
==========================================
  Files         316      327      +11     
  Lines       41948    42603     +655     
  Branches    16853    17143     +290     
==========================================
+ Hits        28498    28978     +480     
- Misses      11186    11345     +159     
- Partials     2264     2280      +16

Impacted Files	Coverage Δ
include/cantera/numerics/FuncEval.h	`57.89% <0.00%> (-26.73%)`	⬇️
include/cantera/thermo/Phase.h	`83.67% <ø> (ø)`
include/cantera/zeroD/Reactor.h	`59.09% <0.00%> (-5.91%)`	⬇️
include/cantera/zeroD/ReactorNet.h	`76.66% <ø> (ø)`
include/cantera/zeroD/MoleReactor.h	`11.11% <11.11%> (ø)`
src/numerics/FuncEval.cpp	`29.16% <14.28%> (-20.84%)`	⬇️
include/cantera/numerics/PreconditionerBase.h	`23.07% <23.07%> (ø)`
include/cantera/numerics/Integrator.h	`15.38% <40.74%> (+10.69%)`	⬆️
src/zeroD/MoleReactor.cpp	`50.00% <50.00%> (ø)`
src/numerics/AdaptivePreconditioner.cpp	`68.33% <68.33%> (ø)`
... and 18 more

📣 Codecov can now indicate which changes are the most critical in Pull Requests. Learn more

ischoegl · 2021-05-27T21:32:08Z

@anthony-walker ... thanks for adding the tests. As a quick heads-up, please avoid doublereal in new code, and keep line lengths to 88 characters (see code style guide / CONTRIBUTING.md). I won't comment on the implementation itself as this is in draft.

speth

@anthony-walker, thanks for the work on this. I think a lot of people are looking forward to this feature. I have some (hopefully helpful) questions about the mathematical / algorithmic issues and suggestions about coding style.

First, as a point of clarification, at the end of AdaptivePreconditioner::setup, it looks like the m_matrix variable is meant to represent the Jacobian. Is this correct? Assuming that it is, I tried printing out the Jacobian for the initial timestep of the Python test problem, and what I saw didn't make a lot of sense. For instance, there are a lot of zeros on the diagonal, entire rows of zeros, and every value is positive. Just as a general sanity check, I would expect the diagonal elements to all be negative (except for a species with no reactions), and the matrix to have a mix of positive and negative values. I tried comparing this Jacobian to a simple finite difference approximation for just the species equation, and saw basically no correspondence.

Also, I don't see anywhere where the preconditioner matrix P = I - gamma * J gets formed. Are you just using the Jacobian directly as the preconditioner? There are places where the two seem to be conflated in the code, but the Jacobian itself presumably isn't what you want to use as the preconditioner.

As a more minor point, I see that the factorization step is being done as part of the "AdaptivePreconditioner::solve" function. To improve performance, the usual strategy is to factorize the preconditioner as part of the "setup" function, so that the repeated calls to "solve" can re-use the LU factors and solving is just a matter of backsubstitution. Later, once this is all working, you might also want to use one of the "incomplete" LU factorization algorithms that Eigen as a way of speeding things up even more, since the preconditioner doesn't need to be exact.

I also had a couple of suggestions on coding style and structure, which will hopefully make the code easier to read, review, and maintain:

There's a lot of manual memory management going on, which I think should be avoided. This means using std::vector instead of new[], std::unique_ptr or std::shared_ptr in other cases where you can't allocate an object on the stack. Basically, you want to avoid ever having to directly use delete.
Inside a member function, it is almost never necessary to prefix access to a member function or variable with this->.
Please take a look at some of the other source code, e.g. Reactor.cpp for typical style and try to write your code to be consistent with this. Take a look at things like indentation (4 spaces, but none for the initial namespace declaration); spacing and placement of brackets in if / else blocks, for loops, etc.; spacing around operators (one on each side of +, = and others).

ischoegl

@anthony-walker ... great to see this moving along. I have one (tangential) comment at the moment.

data/methane_onestep.yaml

ischoegl

Another question.

src/numerics/AdaptivePreconditioner.cpp

include/cantera/numerics/AdaptivePreconditioner.h

ischoegl

Hi @anthony-walker ... thank you for your continued work on this proof of concept! - this is an extremely promising addition that will have tremendous impact once it is fully implemented.

One thing that I wanted to bring up at this point in somewhat more detail is that I believe that the API would benefit from being better integrated with the existing Reactor and ReactorNet implementations (I have hinted at this earlier, see here and here). As you correctly imply, #1089 implemented derivatives directly within the Kinetics object which means that they are fully integrated. Similarly, I believe that the calculation of Jacobians should be fully integrated with zeroD objects and be accessible from their APIs.

Specifically, I think the following (proposed) methods are worth spending some thought on:

virtual Reactor::jacobian (C++) and Reactor.jacobian (Python): while it is ok to just implement these methods for IdealGasMoleReactor / IdealGasConstPressureMoleReactor at the moment, it should be technically feasible to calculate them for other Reactor specializations as well. This would be especially useful if Cantera/enhancements#91 gets tackled (which I believe is a worthwhile effort) and open the preconditioning up to other specializations without having to understand the implementation of the preconditioner (it suffices to implement the correct Jacobian for a given Reactor specialization). It will also make it much easier to build preconditioners based on different versions of Jacobians (i.e. different approximation levels).
ReactorNet::jacobian (C++) and ReactorNet.jacobian (Python): same reasons as above.

I honestly believe that this approach will be tremendously helpful for unit testing: based on a numerical solution, the ODE terms at t + deltat should be close to the ODE term at time t plus Jacobian times difference in states. Imho, this is a much more rigorous test than comparing the results of a time integration.

In addition, this would allow for intuitive interactions with Jacobians, rather than working with the preconditioner directly as currently implemented in the Reactor::preconditionerSetup approach. From my perspective, it would make more sense if the preconditioner would be isolated in ReactorNet and not appear elsewhere.

speth · 2022-04-21T15:59:02Z

@ischoegl, I agree that ultimately, it would be very useful to have a general interface to Reactor and ReactorNet Jacobian matrices. However, this PR already has a lot of things going on, and adds some new features that are immediately useful to users, so I think we should aim to merge it based on the existing feature set, once it's gone through a little bit of cleanup but without necessarily any further major refactoring. Not least of all, it will be much easier to review a future PR refactoring the code to introduce Reactor::jacobian(...) etc. if that's done on its own rather than mixed in with this already large PR.

I'd suggest that for this PR, there are two main things to do to make such a future enhancement as simple as possible:

provide a simple interface for enabling the use of the preconditioner on a reactor network. @anthony-walker has already suggested as net.preconditioner = ct.AdaptivePreconditioner() for the Python interface, and I think we wouldn't need to make much if any change to this.
Mark new methods that introduce coupling between the preconditioner and individual reactors (mainly, this is Reactor::preconditionerSetup) as as "experimental", so we can change or remove them without having to go through a round of deprecations.

anthony-walker · 2022-04-21T16:15:29Z

@speth This is already done.

provide a simple interface for enabling the use of the preconditioner on a reactor network. @anthony-walker has already suggested as net.preconditioner = ct.AdaptivePreconditioner() for the Python interface, and I think we wouldn't need to make much if any change to this.

Is there a specific way to mark methods as experimental?

Mark new methods that introduce coupling between the preconditioner and individual reactors (mainly, this is Reactor::preconditionerSetup) as as "experimental", so we can change or remove them without having to go through a round of deprecations.

ischoegl · 2022-04-21T16:28:04Z

@speth ... I respectfully disagree on the status, but am happy to engage in a discussion (which I've been trying to start for some time).

@anthony-walker ... I honestly believe that some refactoring would make your work a lot stronger and easier to extend/maintain. Is this all the work you can do on this or will there be follow-up PR's?

anthony-walker · 2022-04-21T16:34:04Z

@ischoegl I intend to do follow up PRs, this one has just been drawn out for some time and is becoming massively cluttered and more complicated. I have some quick ideas for refactoring based on your suggestions that may not be too complicated. I will implement them quickly and we can discuss.

@anthony-walker ... I honestly believe that some refactoring would make your work a lot stronger and easier to extend/maintain. Is this all the work you can do on this or will there be follow-up PR's?

ischoegl · 2022-04-21T16:44:30Z

@anthony-walker ...

I intend to do follow up PRs, this one has just been drawn out for some time and is becoming massively cluttered and more complicated.

That is a relief! As long as this is the case, I'm overall not opposed to wrapping this PR up so we can get past this stage. Before doing this, could you post a new work-in-progress issue on Cantera/enhancements (or even in Cantera/enhancements#100)? When I did the refactoring of reaction rate evaluators in Kinetics, I had a master document that listed all the changes I was planning (see Cantera/enhancements#87). In your case, the task should be simpler, but will allow to see the big picture of what you have in mind.

speth · 2022-04-21T17:45:10Z

could you post a new work-in-progress issue on Cantera/enhancements (or even in Cantera/enhancements#100)?

I'd suggest that documenting this plan an existing enhancement makes the most sense. Cantera/enhancements#86 (Provide Jacobian to Sundials solvers) is also very relevant.

Is there a specific way to mark methods as experimental?

Elsewhere, what we've done is to add the following to the Doxygen docstring:

     *  @warning  This method is an experimental part of the %Cantera API and
     *      may be changed or removed without notice.

ischoegl · 2022-04-21T18:07:53Z

@speth ... the original issue is Cantera/enhancements#47 (I believe)

anthony-walker · 2022-04-21T19:45:08Z

@speth @ischoegl It was relatively straight forward to remove the preconditionerSetup methods in lieu of a jacobian method. I also made the preconditioner a part of the Integrator which I think is more intuitive. This allows for just a preconditionerSetup function within the network and not within each reactor specifically.

anthony-walker · 2022-05-06T18:41:18Z

@speth @ischoegl Thank you for the feedback thus far. I believe I have addressed all current comments and concerns and that it is ready for another round of review.

speth

Thank you for the updates so far, @anthony-walker. I have another set of suggestions that I think will continue to simplify things and bring this into line with the other parts of the Cantera codebase. Beyond the various specific comments below, two general requests I have are:

Please run scons doxygen, look at the warnings issued, and fix the ones that are associated with your changes. There are quite a few places where there are discrepancies between the documented parameters and the actual function parameters.
Please go through your code and wrap lines that are longer than 88 characters. I would recommend setting a "ruler" in your IDE to this length so you can see where you are.

include/cantera/numerics/CVodesIntegrator.h

include/cantera/zeroD/IdealGasMoleReactor.h

include/cantera/numerics/FuncEval.h

include/cantera/zeroD/ReactorNet.h

interfaces/cython/cantera/test/test_reactor.py

test/zeroD/test_zeroD.cpp

anthony-walker · 2022-06-07T18:11:23Z

@speth Thank you for the feedback. Hopefully we are getting close to wrapping this up. I believe I have address all requested changes. I set my VS-Code ruler to be 88 characters, I was using Rewrap to wrap lines but that messes with the license statement for some reason. So I went through and manually wrapped any lines over 88 characters.

Please run scons doxygen, look at the warnings issued, and fix the ones that are associated with your changes.

I have also done this and corrected the docstrings.

I did also refactor the applyOptions interface and deprecate the int constants. So you may want to take a look at that.

speth

Thanks for making the previous round of updates, @anthony-walker. I think this is looking pretty good from a structural standpoint. Most (though not all) of the issues I've commented on here are just formatting and documentation concerns.

include/cantera/numerics/Integrator.h

include/cantera/zeroD/ReactorNet.h

src/zeroD/IdealGasConstPressureMoleReactor.cpp

test/zeroD/test_zeroD.cpp

interfaces/cython/cantera/examples/reactors/preconditioned_integration.py

Adds PreconditionerBase and derives AdaptivePreconditioner from this base class. PreconditionerBase is intended to be an abstract class so other types of preconditioners can be extended. AdaptivePreconditioner is a specific type of preconditioning based on a tolerance. The preconditioned integrator works with all Sundials versions supported by Cantera. The ILUT algorithm is used to limit fill-in when factorizing the preconditioner. Adds access to CVODES integrator statistics in C++ and Python Adds reactor classes "IdealGasConstPressureMoleReactor" and "IdealGasMoleReactor" that use a mole-based state vector to increase the sparsity of the preconditioner.

speth · 2022-07-18T00:18:45Z

@anthony-walker - Thanks for the updates. I think this has reached a state where it's ready to be merged, even though there are a few things that need to be done to improve robustness. I just pushed to your branch to squash some of the commits together because there was a lot of churn related to the several times that this PR has been rebased, and the intermediate commits had become impossible to even compile.

@ischoegl - This is awaiting your approval before I can merge it.

anthony-walker · 2022-07-18T03:42:25Z

@speth I agree that it could still use some work to improve robustness. I think said work will be better suited to smaller pull requests now that the core functionality is present.

ischoegl · 2022-07-18T14:15:14Z

@anthony-walker ... could you add what is left to do to Cantera/enhancements#47 (or whatever is most appropriate) so we have a road map for remaining work?

anthony-walker · 2022-07-18T20:53:55Z

@ischoegl yes I will update the most appropriate enhancement.

- Adds a check for surfaces - adds updateState and updatePreconditioner in proper locations - Fixes documentation and makes other formatting changes - Adds utilities tests and expected failures tests for the preconditioner

anthony-walker · 2022-07-19T21:14:59Z

@speth is there anything else that needs done before this is merged?

speth · 2022-07-19T21:48:37Z

@anthony-walker - Just the requested update to Cantera/enhancements#47.

This was referenced Apr 9, 2021

Adding adaptive preconditioner functionality to Cantera (draft) #951

Closed

Refactor Reaction Rate Evaluators Cantera/enhancements#87

Closed

ischoegl added enhancement Performance Reactors labels Apr 9, 2021

bryanwweber marked this pull request as draft May 20, 2021 14:21

bryanwweber changed the title ~~Updated: adding adaptive preconditioner functionality to Cantera (draft)~~ Adding adaptive preconditioner functionality to Cantera May 20, 2021

anthony-walker force-pushed the ap-dev branch 2 times, most recently from 65924db to 2d77873 Compare June 9, 2021 20:10

speth reviewed Jun 11, 2021

View reviewed changes

anthony-walker force-pushed the ap-dev branch from 81656b5 to beb29a5 Compare June 16, 2021 18:32

anthony-walker force-pushed the ap-dev branch 8 times, most recently from ca2294a to 37ccaf5 Compare July 19, 2021 19:31

anthony-walker force-pushed the ap-dev branch 2 times, most recently from 8666869 to 9f28061 Compare July 22, 2021 06:10

anthony-walker force-pushed the ap-dev branch from 834de15 to 64d2b52 Compare July 29, 2021 15:57

ischoegl reviewed Jul 29, 2021

View reviewed changes

data/methane_onestep.yaml Outdated Show resolved Hide resolved

anthony-walker force-pushed the ap-dev branch from 1fb0ed0 to ed4b451 Compare July 30, 2021 14:48

ischoegl reviewed Jul 30, 2021

View reviewed changes

src/numerics/AdaptivePreconditioner.cpp Outdated Show resolved Hide resolved

include/cantera/numerics/AdaptivePreconditioner.h Outdated Show resolved Hide resolved

anthony-walker force-pushed the ap-dev branch from 9039024 to 28cd469 Compare July 30, 2021 17:45

ischoegl requested changes Apr 21, 2022

View reviewed changes

anthony-walker requested a review from ischoegl April 21, 2022 19:45

speth previously requested changes May 29, 2022

View reviewed changes

ischoegl mentioned this pull request Jun 6, 2022

Steady-State Solver for ZeroD Reactors #1021

Closed

4 tasks

speth requested changes Jul 14, 2022

View reviewed changes

speth reviewed Jul 17, 2022

View reviewed changes

interfaces/cython/cantera/examples/reactors/preconditioned_integration.py Outdated Show resolved Hide resolved

anthony-walker added 2 commits July 17, 2022 19:38

Change temperature tolerance in integrate

d232f6d

speth previously approved these changes Jul 18, 2022

View reviewed changes

ischoegl previously approved these changes Jul 18, 2022

View reviewed changes

Address PR review comments

e74ffd3

- Adds a check for surfaces - adds updateState and updatePreconditioner in proper locations - Fixes documentation and makes other formatting changes - Adds utilities tests and expected failures tests for the preconditioner

speth approved these changes Jul 19, 2022

View reviewed changes

anthony-walker mentioned this pull request Jul 20, 2022

Addition of Jacobian Preconditioning Functionality to Cantera. Cantera/enhancements#47

Open

speth mentioned this pull request Aug 26, 2022

Preconditioned solver improvements #1377

Merged

5 tasks

Uh oh!

Conversation

anthony-walker commented Apr 9, 2021 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

anthony-walker commented Apr 9, 2021

Uh oh!

ischoegl commented Apr 9, 2021

Uh oh!

codecov bot commented Apr 9, 2021 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Codecov Report

Uh oh!

ischoegl commented May 27, 2021

Uh oh!

speth left a comment

Choose a reason for hiding this comment

Uh oh!

ischoegl left a comment

Choose a reason for hiding this comment

Uh oh!

Uh oh!

ischoegl left a comment

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

ischoegl left a comment • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

speth commented Apr 21, 2022

Uh oh!

anthony-walker commented Apr 21, 2022 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

ischoegl commented Apr 21, 2022

Uh oh!

anthony-walker commented Apr 21, 2022 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

ischoegl commented Apr 21, 2022 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

speth commented Apr 21, 2022

Uh oh!

ischoegl commented Apr 21, 2022

Uh oh!

anthony-walker commented Apr 21, 2022

Uh oh!

anthony-walker commented May 6, 2022

Uh oh!

speth left a comment

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

anthony-walker commented Jun 7, 2022 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

speth left a comment

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

anthony-walker commented Apr 9, 2021 •

edited

Loading

codecov bot commented Apr 9, 2021 •

edited

Loading

ischoegl left a comment •

edited

Loading

anthony-walker commented Apr 21, 2022 •

edited

Loading

anthony-walker commented Apr 21, 2022 •

edited

Loading

ischoegl commented Apr 21, 2022 •

edited

Loading

anthony-walker commented Jun 7, 2022 •

edited

Loading