StoichManager vs Eigen

[ischoegl]

Looking ahead to some work that we're planning in the near future, I was wondering about second opinions for replacing StoichManager by sparse matrices. Presumably, Eigen should be capable of handling most operations out-of-the-box, and would allow for a seamless integration with std::vector, among other benefits (it would simplify adding&removing of species and reactions, but there are other things I can think of as well). In terms of speed, some tests would have to be run, but any optimization that is part of Eigen would be automatically leveraged. Edit: probing a little deeper, it appears that most of the operations I had looked at are only available for dense matrices; on the other hand, you can still build on the library; so this may be more involved.

Any comments/thoughts would be appreciated (@speth, @kyleniemeyer, et al).

[bryanwweber]

One of the things I think we can consider here is developer ergonomics. On one hand, sparse arrays may be more intuitive for those folks familiar with linear algebra methods; on the other hand, the specialized StoichManager may be more familiar for those from a chemistry background. I suspect that using standard array/matrix methods for calculations (applied to sparse matrices) would be easier for most developers (myself being one), but I'd need to see a bit of the implementation to really think about it deeply. Sorry, this is kind of a non-answer 😂

[ischoegl]

See Cantera/cantera#1081 ... ended up testing a hybrid solution (keeping StoichManager for the gnarly stuff, while outsourcing standard sparse matrix multiplications to Eigen). Code is much shorter / easier to read ... 🎉

[ischoegl]

It's slightly slower though. Still, this led to some stuff we have been interested in for a while (sparse Jacobians) ...

[speth]

I think using Eigen's SparseMatrix class as the container for working with and providing access to the Jacobian makes sense. But as I think you've already seen, a lot of the low level operations that need to be done don't quite correspond to operations that need to be done by the StoichManagerN, which takes care of some tricky corner cases in how to calculate rates when given negative species concentrations and has a lot of optimizations based on the fact that the non-zero elements of the stoichiometric coefficient matrix are mostly ones with the occasional two thrown in.

In what you currently have in Cantera/cantera#1081, I don't quite understand what using Eigen is accomplishing outside of the functions used in calculating the Jacobian. The new dotReaction and speciesDot methods both introduce additional copying steps every time they are called, and which were unnecessary in the previous code.

[ischoegl]

But as I think you've already seen, a lot of the low level operations that need to be done don't quite correspond ...

Absolutely!

I agree that both have advantages and StoichManagerN is certainly not going away. Buy as you mention, sparse matrices are useful for the calculation of Jacobians, as forward and reverse directions aren’t aware of each other’s coefficients, and the multiplication is much easier. Regarding dotReaction and speciesDot - those don't copy in the newer versions in Cantera/cantera#1081 as some of those are just writing to MappedVectors (at least that’s the idea), but there is one more addition while avoiding zeroing out the buffer initially. The old implementation keeps aggregating/accumulating, which is highly optimized (and faster than Eigen), but doesn’t lend itself to a calculation of partial results. Having almost duplicate code for addition/subtraction should be avoidable if MappedVectors are used, but it's not something I'd fight over.

What I am trying to do at the moment is to keep clean versions of individual terms in the product so they don’t have to be recalculated for the Jacobian. Eigen’s algebraic capabilities are quite neat here (loops are largely avoided) and should make for easy-to follow code. (I am currently staring at calcROP which is where all the aggregation takes place.) Once I have this working, I can bench-mark.

On a separate note, I have already started talking to @anthony-walker … . It appears that the new work blends nicely with Cantera/cantera#1010 (also, first speed comparisons are promising).