Add Power Electronics sparsity pattern computation into allocate

[alexander-novo]

Description

Adds sparsity pattern computation into PowerElectronicsModel::allocate().

Closes #301

Checklist

Put an x in the boxes that apply. You can also fill these out after creating
the PR. If you're unsure about any of them, don't hesitate to ask. We're here
to help! This is simply a reminder of what we are going to look for before
merging your code.

• All tests pass.
• Code compiles cleanly with flags -Wall -Wpedantic -Wconversion -Wextra.
• The new code follows GridKit™ style guidelines.
• There are unit tests for the new code.
• The new code is documented.
• The feature branch is rebased with respect to the target branch.
• I have updated CHANGELOG.md to reflect the changes in this PR. If this is a minor PR that is part of a larger fix already included in the file, state so.

Further comments

I added a couple of helper structs, but these are only there for cache locality and readability purposes. They are local to the function itself. Additionally, I made use of std::vector as a dynamically growing container, but all final products end up as simple pointers.

[alexander-novo]

I also wanted to add an additional comment on the removal of the short-circuiting condition in COO_Matrix::axpy():

I found out yesterday that IDA sets the alpha parameter to 0 during IDACalcIC(), which was causing the sparsity pattern to change. This (for some unknown reason) was what was causing KLU to segmentation fault. So this change also fixes the issue that @abdourahmanbarry noticed in his fault simulation.

[pelesh]

Great progress. Some design tweaks are needed but we are almost there..

[pelesh]

I would avoid using STL containers here. What we need to do is to

• populate row pointers and column indices in the CSR Jacobian matrix
• have components store mapping from their Jacobian elements to the index of column indices array

[alexander-novo]

I cannot really remove the STL container for reverse_map... the size of the second level arrays are unknown ahead of time and need to be dynamically grown. I can remove the vector for component_sparsities, though

[pelesh]

I am not sure this is a good idea. Ideally, one should loop through components and components would read/write at locations determined by the connectivity information. The system composer should remain agnostic of specific details of the components.

[alexander-novo]

Well a component cannot write all of its info at once, so I must have some way of differentiating the different parts of info that a component needs to write at once.

[alexander-novo]

rebase w.r.t. develop

[alexander-novo]

I changed the sparsity to be calculated into CsrMatrix object, but because of the aforementioned issue with the assignment operators' implementations in CsrMatrix, there is a use-after-free bug.

@nkoukpaizan maybe it's worth fixing? I can do some wonky stuff to try and avoid it but it seems to me like maybe it's better to just fix the bug?

[nkoukpaizan]

I changed the sparsity to be calculated into CsrMatrix object, but because of the aforementioned issue with the assignment operators' implementations in CsrMatrix, there is a use-after-free bug.

@nkoukpaizan maybe it's worth fixing? I can do some wonky stuff to try and avoid it but it seems to me like maybe it's better to just fix the bug?

Wouldn't it work if we stored csr_jac_ as a pointer, using new and deallocating in the system model's destructor? My understanding is that we don't want to use the assignment operator.

[alexander-novo]

I changed the sparsity to be calculated into CsrMatrix object, but because of the aforementioned issue with the assignment operators' implementations in CsrMatrix, there is a use-after-free bug.
@nkoukpaizan maybe it's worth fixing? I can do some wonky stuff to try and avoid it but it seems to me like maybe it's better to just fix the bug?

Wouldn't it work if we stored csr_jac_ as a pointer, using new and deallocating in the system model's destructor? My understanding is that we don't want to use the assignment operator.

That's the wonky stuff I meant. Is CsrMatrix only ever meant to exist on the heap, then? That should probably be mentioned somewhere.

[pelesh]

Wouldn't it work if we stored csr_jac_ as a pointer, using new and deallocating in the system model's destructor? My understanding is that we don't want to use the assignment operator.

I don't see are need to use assignment operator either (and it might be a good idea to delete it tbh). The implementation should go something like this:

1. In constructor csr_jac_ is instantiated.
2. In allocate function, NNZ is computed. After that, the Jacobian is allocated with the computed NNZ.
3. Next (also in allocate), the sparsity pattern of the Jacobian is computed.

This is what is needed to address #301.

[alexander-novo]

Wouldn't it work if we stored csr_jac_ as a pointer, using new and deallocating in the system model's destructor? My understanding is that we don't want to use the assignment operator.

I don't see are need to use assignment operator either (and it might be a good idea to delete it tbh). The implementation should go something like this:

1. In constructor `css_jac_` is instantiated.

2. In `allocate` function, NNZ is computed. After that, the Jacobian is allocated with the computed NNZ.

3. Next (also in `allocate`), the sparsity pattern of the Jacobian is computed.

This is what is needed to address #301.

But as far as I can tell, there is no way to change the size of the matrix after construction. And the size of the system is not known at construction time for PowerElectronicsModel

[pelesh]

That's the wonky stuff I meant. Is CsrMatrix only ever meant to exist on the heap, then? That should probably be mentioned somewhere.

You could also have

class SystemSolverPowerElectronics
{
public:
  // some code

private:
  CsrMatrix jacobian_;
};

In that case Jacobian matrix would be instantiated with a default constructor and you would allocate it and populate its data arrays. That wouldn't change much in the approach.

[pelesh]

But as far as I can tell, there is no way to change the size of the matrix after construction. And the size of the system is not known at construction time for PowerElectronicsModel

You implemented CsrMatrix class, so feel free to add methods that may be needed for the system Jacobian assembly. 😄 Adding a resize method would make sense to me in this case.

[alexander-novo]

I copied the implementation from ReSolve, I'm not really familiar with how it works. It kind of seemed to me like there may have been a good reason why there wasn't a resize() method. Would a resize() method that simply sets m,n be correct?

[alexander-novo]

Well, as far as I can tell m_ is used for nothing and n_ is just used for allocating memory, so I think it should probably be harmless to just modify them. Please let me know.

[alexander-novo]

Okay, I went ahead an imposed the correct system vector ordering. The Microgrid examples already have this ordering, but I did have to change the RLCircuit example to have the correct ordering.

Now that the ordering is required to be nice, I went ahead and updated the jacobian sparsity computations to take advantage of it - first the internal rows are done component-by-component, then the heavier reverse-mapping machinery is invoked only for the external rows.

I also realized that I could use the sum of the component Jacobian nnzs as an upper bound for the system Jacobian nnz value and got rid of the vector for global_col_indices by just pre-allocating this upper bound. I don't think that this will be too big an over-estimate to be concerned about the wasted memory.

[alexander-novo]

@pelesh I implemented what we discussed at the SciDAC meeting, and am not depending on the components having their jacobians in CSR form anymore (just using COO).

The two outstanding pieces of feedback, though, I'm waiting for your feedback on. Maybe we can talk about them in a code review?