Hi, I'm currently playing with implementing (proper) support for petsc in OPM. For this we need better containers and interfaces than petsc's own.
The thing is, petsc doesn't support direct memory access. While this is good object oriented practice, it creates some headaches when implementing references to the scalar in an arbitrary cell. The question is: do we want at()/operator[][] for matrices at all? Or is it sufficient to say "hey pal, if you want to access the value in a specific cell you need to go through iterators"?. While it is possible to implement reference-esque features towards petsc, it will be subject to data races. I looked up deal.II, a FEM library that can use petsc, and they do not support random access for matrices. http://www.dealii.org/8.2.0/doxygen/deal.II/classPETScWrappers_1_1BlockSparseMatrix.html Regarding data races: Iterators are needed, that's a no-brainer. I think of doing this the deal.II way, i.e. reading per-row into a buffer which in turn writes to the matrix in the case of operator*. This again can either lead to a data race if there is no global look-up mechanism for uniquely spawned instances to share buffers (which again provides races :)), or we can simply declare it undefined behaviour to have two live iterators to the same row. Since petsc uses much MPI this isn't too unreasonable. Of course, we can extend the same "this-is-undefined" statement to also concern scalar references, but I personally think that this much undefined behaviour of unrelated means of access is just bad design, and I'd rather omit it. So: do we need or want per-cell references to matrices? Finally, some code examples. Keep in mind that petsc ONLY supports GetRow(). Consider this case of iterators. Perfectly defined and as you expect. matrix m; for( auto itr = m.begin(); itr != m.end(); ++itr ) *itr += 3; Now, another well-formed example. Notice that other_itr was created through itr. They share underlying buffers. for( auto itr = m.begin(); itr != m.end(); ++itr ) { auto other_itr = itr; while( other_itr++ != m.end() ) { *other_itr += 3; } } An undefined example: auto itr1 = m.begin(); *itr1 += 3; auto itr2 = m.begin() + 2; *itr2 += 4; ++itr1; *itr1 += 2; // POSSIBLY UNDEFINED. *itr1 += 2 may or may not be undefined (by us). It's ok in the sense that the underlying matrix isn't written in this cell by itr2, which is ok. However, assuming writes aren't deferred itr2's buffers up until line 4 are identical to itr1's. Should itr1 increment once more it will not have received writes from itr2 and reading would definitely be undefined. This can be partially solved with subscription and (global) cache invalidation, but personally I think this adds more complexity than gain, so I'm with non-enforced read-on-write being undefined. Unfortunately petsc introduces the issues of concurrent read and write even in the single-CPU single-threaded case, and provides no elegant means of preventing it.
signature.asc
Description: OpenPGP digital signature
_______________________________________________ Opm mailing list [email protected] http://www.opm-project.org/mailman/listinfo/opm
