Hi Jørgen, <disclaimer> I am one of the guys behind DUNE and especially dune-istl. Therefore I might have wrong stereotypes concerning other packages and be not objective in this regard. </disclaimer>
On Tue, Oct 28, 2014 at 12:14:47PM +0100, Jørgen Kvalsvik wrote: > 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. Thanks for nice introduction. I like that as I prefer to know who I am dealing with. > 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. > I admit that I am not exactly aware what the upscaler is. But please be aware that there also is another LinearSolverInterface in opm-porsol. > 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. I beg to differ here (and agree with Arne). One should not try to see linear algebra as sparse matrices and vectors with some added sugar. Actually it is the sugar (read the solvers) that is important. In the real world you either decide to use a package (which you probably already know) or more ideally you look for the best kind of solver, choose a package, and use its data structures. If you want flexibility then you have to take into account that each package might use its own data structure (mostly for good reasons, sometimes not). Even if they use the same format, setting up the data structure might be quite different (e.g. PetSc vs. dune-istl). With the rise of new hardware architectures the sparse matrix format might even be different depending on the platform, e.g. coalesced CSR or Ellpack for GPUs is quite common. If you limit yourself to one package per simulation run, this means that users must be able to choose the package and the algorithms need to adapt to this. In dune-pdelab <https://github.com/dune-project/dune-pdelab/tree/master/dune/pdelab/backend> (not used in OPM) we have different backends for (sequential) PetSc, Eigen, and dune-istl that provide a common way of setting up matrices and vectors during the discretization and provides a solver interface. Basically each solver is a class and has an apply method, see e.g. <https://github.com/dune-project/dune-pdelab/blob/master/dune/pdelab/backend/eigen/solvers.hh>. Just supporting these 3 packages is a huge effort, which might not even payoff and changes to the codebase often break some packages (apart from dune-istl of course :). If you are targeting parallel solvers, the picture gets even more complicated as each package of course has its own way of dealing with parallelism in datastructures :( > 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. That is indeed a downside we should eliminate. Then again we want to rule the world ;) > > 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. > Like I said, I highly doubt that such a format exists (at least when it comes to parallel solvers) or that we can persuade other projects to use 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. > Using C might for a reason. There several parts (at least within opm-core) that must/should be usable from C-code if I am not mistaken. > 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). > The *mess* might also be for a reason. Namely, efficiently setting up the data structures and filling them. Your are pushing triples (row,col,val) into a vector of unknown size. Won't this result in reallocating memory and copying old data whenever the vector realizes that it is too small? Maybe having the possibility to provide an upper bound of nonzeros could help. Ease of use also comes for price when calling to_csr which needs quite some passes, data movements... Please be sure to compare the performance here. Some other things/tips that might be worth investigating: - Maybe you should make the row_ind_ array one entry larger to prevent branching when accessing elements. I think this is standard for most CSR implementations - Run guess_cols only once. - Use std::lower_bound instead of handed coded binary search on sorted arrays. Cheers, Markus -- Dr. Markus Blatt - HPC-Simulation-Software & Services http://www.dr-blatt.de Hans-Bunte-Str. 8-10, 69123 Heidelberg, Germany, USt-Id: DE279960836 Tel.: +49 (0) 160 97590858
signature.asc
Description: Digital signature
_______________________________________________ Opm mailing list [email protected] http://www.opm-project.org/mailman/listinfo/opm
