28. okt. 2014 kl. 12:14 skrev Jørgen Kvalsvik <[email protected]>:
> 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. Thank you for introducing yourself and your project, it is always good to know what people are doing with OPM! > > 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. This is true. I have considered replacing the direct use of Dune::BCRSMatrix in the upscaler's flow solver (class IncompFlowSolverHybrid) with using opm/core/linalg/LinearSolverInterface, I have not done so, mostly for the following reasons: - Why change working code? - No resources (time) to do it. - Potential downsides (overhead from matrix conversion etc.) Note that there are more cases of matrix conversion that you may not be aware of. For example, the automatic differentiation code uses Eigen and its SparseMatrix class, in fact a somewhat ugly solution was implemented in opm-autodiff/opm/autodiff/DuneMatrix.hpp to create a Dune compatible matrix from an Eigen matrix. > > 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. I am a bit unsure about what we can really do. Ideally, if we choose (or make) a single such representation, it should have the following properties: 1. Can be built in many ways, to be fit-for-purpose for different types of assembly (FEM, FV, Mimetic etc.). 2. Has minimal memory overhead. 3. Representations suitable for any third-party linear solver we'd like to run can be easily (and quickly) extracted. Combining the first and second point is hard enough, but I think that the third point is the hardest. For example: - UMFPACK requires a CSC, not CSR matrix. - Building a Dune matrix (at least before Dune 2.3) takes a lot more time than just copying the data. - Distributed parallel libraries may require the data to be stored in particular distributed vector and matrix objects. There exists an alternative to having one, common format/class that still gives a lot of flexibiltity: a templated approach, where you directly build the matrix type of your choice (Petsc, Dune etc.) with functions that wrap these with a unified interface. To some degree, we have already tried this (see opm/core/transport/implicit/JacobianSystem.hpp and other files in the same directory). A disadvantage is that the code quickly becomes complex and hard to understand for anyone but experts. Also, we would like to have run-time flexibility, not just compile-time flexibility. > > 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 do not think that the call_*.c files need to be removed, but I do not think new solver interfaces should be required to provide a C interface. I would certainly consider such a proposed SparseMatrix class. It should be placed in opm/core/linalg I think. However, I am afraid that it may end up as "Yet Another Sparse Matrix" and only add to any complexity or confusion there already is. Abandoning the Eigen SparseMatrix for opm-autodiff will require a substantial effort, that I do not think we'll make (at this point at least), so I think that will remain. If your proposal gives us greater flexibility for the upscaling code without hurting performance, I am inclined to accept it anyway, even if that is the only use case for the new class. As a side note, keep in mind that for the moment, we require all code to build with g++ 4.4, which places some limits on C++11 features you can use (for example auto is ok, but lambdas are not). > > 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. I hope that I did not dishearten you! It is an important question, and a hard one. Atgeirr _______________________________________________ Opm mailing list [email protected] http://www.opm-project.org/mailman/listinfo/opm
