On 10/28/2014 08:56 PM, Markus Blatt wrote: > 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> Duly noted.
>> 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. I'll have a look at it as soon as possible. >> 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. Unifying these structures (and preferably make the calling code oblivious to the fact that they're different) could be a nice solution and definitely something to work on. > > 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. It was considering this I started the templated solution the way I did, because I'm still not certain why they should have very different operations on them. Of course, using template specification, adding individual features from the underlying class is simple. > 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 :( It's the common way of setting up matrices I'm trying to provide. It's work to maintain, surely, but that's an unavoidable price to pay for multiple backend support. A benefit is that this maintenance now is focused only on the data structures and solvers, not scattered around the codebase. >> 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 ;) It's perfectly fine as long as the only solver is Dune. >> >> 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. The primary goal is a nicer interface for internal use that works with more backends. > >> 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. > Maybe, but it seems inconsistent with the (call for) removal of call_petsc.c. The calling code itself is mostly used from within C++ again, but of course I see the value of a C-compatible interface. > >> 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. Performance will of course be measured once things starts working. There is a cost to reallocating (and copying) memory during resizes. Typically this is quite small, but supporting a guess-size or specific realloc, basically a forward to std::vector::resize is perfectly possible. Using the triples vector was just an immediate and simple solution. The builder class was also made with different builder backends in mind. Building directly into some compressed format is usually quite expensive, and from I can read in IncompFlowSolverHybrid it a lot of the code is actually resizing of internal buffers in Dune::BCSRMatrix. Using a separate builder completely removes this programmer cost. Of course, I have yet to actually measure the difference because it isn't working yet. When I get a working example I'll perform some testing. > > 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. > Thanks, I'll have a look. I suspect the implementation will be somewhat different when it is completed. Sincerely, Jørgen
signature.asc
Description: OpenPGP digital signature
_______________________________________________ Opm mailing list [email protected] http://www.opm-project.org/mailman/listinfo/opm
