Adding adaptive preconditioner functionality to Cantera

[anthony-walker]

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]

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

[ischoegl]

Tagging #951.

[codecov]

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]

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

[ischoegl]

Another question.

dismissing to allow re-request

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

[speth]

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

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

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

@ischoegl yes I will update the most appropriate enhancement.

[anthony-walker]

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

[speth]

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