Petsc interface generalisation

[cmacmackin]

As described in #1873, this generalises the PETSc interface so it can be used on a wider variety of problems. This PR contains two significant changes:

• the introduction of an OperatorStencil class which provides information on a matrix's sparsity pattern, needed for preallocating memory
• the OperatorStencil is further leveraged to determine which guard cells should be given unique global indices and included in PETSc vec objects.

These changes will eventually be merged as part of #1803, but as work is ongoing to apply the PETSc interface to other parts of the code I thought it would be best to get them merged into next ASAP.

[ZedThree]

Looks good! I like the new way of getting the "thin" regions.

More methods on the PetscVector/Matrix could be made const. I'd try to make all of them const at first and them remove const until it works.

Please you could apply clang-format to these changes too?

[cmacmackin]

Thanks for the comments. Some of these are duplicates of those on #1803, which I am in the process of addressing right now and will place on this branch once ready.

[cmacmackin]

@ZedThree I'm not sure if many more methods on the vector and matrix classes should be made const, actually. While strictly speaking they can be (as it is not the contents of the class itself which is modified) they still are affecting changes to the values stored in the backend PETSc object. Thus while the compiler wouldn't complain they are semantically non-const.

A question on style/convention when it comes to these issues: is it every appropriate to apply const to the function but not to the function result? Specifically, it would be possible to make the indexing methods const as they only return an accessor for an element of the vector/matrix. However, it makes no sense to declare these accessors as const, as their main purpose is for setting the value of elements. Given that the indexing would usually be used to return an accessor which would immedoately be used to set the value of an element, it strikes me that it would be misleading to declare these methods as const, even though strictly speaking they are. What would normally be done in this case?

[ZedThree]

You're right that's there's not much more on the vector/matrix classes that could be const, I was thinking of GlobalIndexer, apologies! Things like getGlobal and isLocal should be const. Less certain about e.g. getNumDiagonal, as it would need numDiagonal to be mutable. Probably ok if you view the latter as a cache variable. It should probably at least return const& though to avoid a copy.

Indeed, a const method means "logically const" as opposed to "physically const", so it's not really correct to make something const that returns a object that can be used to modify the state of the parent. When it comes to these things, we can add a const overload that returns a const& or const*. I think you already have this for PetscVector::get. Dropping the const from the return type only makes sense for simple types like double or int (really, value types) where they can do harm, as it were.

Element is really a reference type, and as you point out, should be const. However, for the indexing, I had the following overloads for the HypreMatrix::operator():

Element operator()(const ind_type& row, const ind_type& column);
BoutReal operator()(const ind_type& row, const ind_type& column) const;

where the latter just gets the value directly. You could return a const Element, but as the only thing you (probably?) want to do with it is to cast it to a BoutReal, I just skipped that part. Also there was some trickiness to creating a const Element to begin with for reasons. Also also, note the const on the arguments!

I think PetscVector::toField should be const and return just T. Other than that, you're right, there's not much else on vector/matrix to make const.

[cmacmackin]

Yes, unfortunately getNum[Off]Diagonal() can not be const because it caches the result. I have switched it to returning a const referance however, and am not sure how I missed that earlier.

Currently I have methods to convert from an Element to a BoutReal but it probably makes more sense to return the latter directly from operator(). It also allows for the const overload to be simpler than the normal indexing method (at least in the case of matrices), as it doesn't need to worry at all about getting positions and weights.

Yes, I'd already spotted that toField() should be made const.

[cmacmackin]

I believe the latest push should address all of your comments. I think it should pass the CI too, although we'll find out in an hour or so.

[ZedThree]

Just to be clear on the getNum[Off]Diagonal(): it can be made const if num[Off]Diagonal is made mutable. This allows the member to be changed even in const methods: remember, this is ok as long as the method doesn't logically change the instance, which is the case here because num[Off]Diagonal is just acting as a cache. Calling getNum[Off]Diagonal doesn't result in any observable changes to GlobalIndexer*.

Also toField returns const T, when it could be T (const on returned value types is most likely ignored by compiler anyway).

Neither of these are blocking changes and can be fixed later, I just wanted to clarify my earlier points!

* except via a timing side-channel, and we don't care about that here.

[cmacmackin]

There are some issues with the OpenMP builds which I am currently addressing. More concerning are the failures with the clang build, which I am at a loss to explain.

[ZedThree]

The clang build has some warnings about inconsistent or missing virtual/override -- it's worth addressing these, even if they don't fix the problem.

OpenMP I guess you've got sorted, but probably due to the gtest ASSERT_EQ inside BOUT_FOR -- just change to BOUT_FOR_SERIAL

[cmacmackin]

I've fixed those issues but the errors with the clang run still persist. What is even stranger is that one of the OpenMP runs still breaks: the one built with CMake. The autoconf build works just fine, despite apparently using the same build-options as the CMake one. I have no idea what could be causing that problem.

[cmacmackin]

The problem in the clang build proved to be quite subtle. Although I was caching the sparsity pattern calculated in the GlobalIndexer class, I actually recalculated it any time either GlobalIndexer::getNumDiagonal() or GlobalIndexer::getNumOffDiagonal() were called. This could invalidate pointers to data from any previous calls. While that didn't cause problems with GCC, it did with clang. I've now fixed the caching so cached values are actually used when available, without recalculating. That fixes the problem.

[cmacmackin]

I don't fully understand what the issue was with the CMake/OpenMP case. However, it seems that IndexerTest/TestGetRegionBndry was not properly cleaning up after itself, resulting in PETSc not being properly initialised in future tests. The problem seems to have been related to a for-loop with the iterator being a reference to GlobalIndexer objects. Switching to using pointers got rid of the problem (when running locally, at least).

[cmacmackin]

@ZedThree
All tests are passing and all requested changes have been made. I believe it is now ready to be merged.

[ZedThree]

Thanks @cmacmackin !

I suspect the issue with that test was something to do with making a copy of the GlobalIndexers, but I can't quite see what the actual problem is. The fix is fine though, so I won't worry about it too much.