Hi,
I am a student at NTNU http://www.ntnu.no/, and my pre-master project is
modularizing OPM by implementing (at least trying to) support for using
more linear solvers. There is some code available, see the
LinearSolverInterface and friends in Core, but there is very little code
that actually uses this. The biggest sinner I've found so far is the
Upscaler which has dune hard-coded as its solver.
Which brings me to some background on the issue I want to discuss.
Support for more linear solver libraries means we need some control of
the data representation. The upscaler uses Dune's data structures both
for building the problem as well as solving it, but these to me seems
unsuited as they're Dune-specific and hard to convert to some more
suited representation.
There is an ongoing effort to support petsc within OPM (see
https://github.com/qilicun/opm-core/pull/2). In Liu Ming's original pull
request https://github.com/OPM/opm-core/pull/590 the question is raised
why he uses C and why the C++ source file is needed at all. It seems to
me that the same issue actually applies to umfpack too, through its
opm/core/linalg/call_umfpack.c and related sys_sparse.c file.
The sys_sparse.c source file actually introduces the CSRMatrix (C)
struct, a tiny implementation of a CSR sparse matrix data structure. The
LinearSolverInterface actually even accepts this structure, unpacks it
and hands it to the other virtual solve method.
The reason I bring this up is because I think we're at a crossroads
here: either we must select a base matrix representation suitable for
all solvers (or easily convertible to them) from some external library,
or we must write our own.
Obviously, the first solution is easier in the sense that development
and maintenance is offloaded to some external project. However, this
might also mean less flexibility in some areas (such as conversions) and
another dependency for the project. Considering the fact that the CSR
structure itself is quite simple combined with the fact that we already
have a CSRMatrix implementation leads me to think that someone else have
thought about this issue before.
There are issues with the (current) struct CSRMatrix implementation too.
First of all, it is written in C. While this in itself isn't necessarily
a bad thing, its interface feels clumpy and difficult and doesn't really
align well with the rest of the project's style. In addition, in light
of the removal of call_petsc.c, I also want to remove call_umfpack.c and
replace it with a nicer C++ implementation. This removes the current
usage of struct CSRMatrix.
The last week I've been working on a C++ SparseMatrix implementation.
It's written in pure (templated) C++. To make the interface nicer I've
also split the matrix building process into its own class, which solves
the state mess that is prevalent both in struct CSRMatrix and the Dune
*CSR matrix classes. Under the hood it is implemented as plain C arrays,
but with the RAII idiom in place. You can see the WIP in my csr-matrix
branch https://github.com/jorgekva/opm-core/tree/csr-matrix with an
emphasis on WIP. I've constructed the class so that other backends
easily can be written and used, both in construction (e.g. list of
lists, dictionaries etc) as well as representation (BCSR comes to mind).
What I propose is that the call_*.c files are removed in favour of C++
implementations. I also propose that my SparseMatrix gets a place
alongside SparseVector and SparseTable in core/utiliy so we get a
uniform and flexible CSR implementation we can program for that is
easily convertible to whatever format the linear solvers need.
I am in the process of porting the IncompFlowSolverHybrid from Dune
matrices to my own format to try this out. I would, however, like some
feedback on this work and if it is the right way to go or if the there
are some other ideas on how to solve this issue.
Sorry for the wall of text, but I felt some background was necessary.
I look forward to the discussion.
Sincerely,
Jørgen Kvalsvik
_______________________________________________
Opm mailing list
[email protected]
http://www.opm-project.org/mailman/listinfo/opm
