Bug: Configuring sparse Jacobians with IDA assumes Jacobian has been evaluated already

[alexander-novo]

Description how to reproduce the bug

Design a simulation which uses sparse Jacobians and does not evaluate the Jacobian before running the simulation (for an example of such a simulation, see this Microgrid benchmark). The Jacobian will not be calculated correctly, and the linear solver in SUNDIALS will error, complaining about being configured incorrectly.

GridKit™ version

develop

System and environment details

N/A

Additional information

The reason for the bug is simple:

This commit changed the sparse linear solver configuration code in ida.cpp to no longer assume that there were n^2 non-zero elements in the jacobian. To do that, it calls (model_->getJacobian()).nnz() but getJacobian() must be called only after calling evaluateJacobian() (un-documented, but true).

The reason that this wasn't caught is because currently, the only simulations we have being tested using sparse Jacobians coincidentally all evaluate their Jacobian before running the simulation:

• See here
• And here
• And also here

But...

1. This shouldn't be a necessary step.
2. Evaluating Jacobians before configuring the simulation isn't even enough - some (potentially) nonzero elements are squashed during Jacobian evaluation if the alpha coefficient is set to 0, which it often is if updateTime() hasn't also been called beforehand.
3. The number of nonzero elements should probably be allowed to change during simulations - discrete events like faults can introduce new variables into the system.

Luckily, SUNDIALS has an interface for re-allocating sparse matrices, which we can call after evaluating the system Jacobian and realising there isn't enough space.

So the correct steps should probably be something along the lines of:

1. When configuring the linear solver, allocate a small (maybe 0 sized if SUNDIALS lets us) jacobian matrix, and keep track of how much space it has for nonzero elements
2. After evaluating a system jacobian, compare the number of nonzero elements it has and how much space the SUNDIALS matrix has - if they differ, then re-allocate the SUNDIALS matrix to fit
3. Finally, copy the entries from the system Jacobian to the SUNDIALS matrix.

Eventually, the GridKit CSR matrix interface should also be able to handle this re-allocation automatically, without having to copy between two buffers.

[alexander-novo]

I implemented a proof of concept fix in my WIP branch:

340384f

It should be easy to add the to develop

[pelesh]

A few comments:

• Constant alpha should not be set to zero. Having alpha == 0.0 means the step size is infinite.
• Jacobian sparsity pattern should not change during the simulation. If the model anticipates discrete event that would add nonzeros, then these should be modeled as structural nonzero and set to zero values.
• The matrix object should have defined nnz_size and nnz_capacity where nnz_capacity > nnz_size, similar as in the case of std::vector. This should allow us to change the nnz without reallocating in an unlikely case when sparsity pattern changes.

[alexander-novo]

• But if you haven't set alpha to anything by calling updateTime() many operating systems will serve you zero-ed out memory, so alpha will just be default 0.
• Sure that's fine, but this is leaving performance on the table, I think.
• Our matrix objects do this, but SUNDIALS as far as I know does not. @Steven-Roberts maybe you can comment on this. I think nnz_capacity must equal nnz_size from what I'm reading in documentation.

[Steven-Roberts]

Our matrix objects do this, but SUNDIALS as far as I know does not. @Steven-Roberts maybe you can comment on this. I think nnz_capacity must equal nnz_size from what I'm reading in documentation.

They can be different. SM_NNZ_S(A) gives the capacity and (SM_INDEXPTRS_S(A))[SM_NP_S(A)] gives how much of that is actually used. SUNSparseMatrix_Realloc(A) will make them equal and SUNSparseMatrix_Reallocate(A, nnz) can make them unequal.

[alexander-novo]

Right, okay, so in SUNDIALS terminology nnz just ends up being capacity. Makes sense.