#854 Improve performance of MCMGMapper
Metadata
| Property | Value |
|---|---|
| Reported by | Martin Nolte (nolte@mathematik.uni-freiburg.de) |
| Reported at | Dec 13, 2010 12:13 |
| Type | Feature Request |
| Version | Git (pre2.4) [autotools] |
| Operating System | Unspecified / All |
| Last edited by | Christoph Grüninger (gruenich@iws.uni-stuttgart.de) |
| Last edited at | Nov 5, 2013 20:58 |
| Closed by | Christoph Grüninger (gruenich@iws.uni-stuttgart.de) |
| Closed at | Nov 5, 2013 20:58 |
| Closed in version | Unknown |
| Resolution | Implemented |
| Comment | Merged to master in …ca709f81aca. Thanks to Martin and Steffen for their ideas, support, and help. |
Description
The finite volume example from the dune-grid-howto can be sped up by more than 5% using the attached version of the MCMGMapper. This implementation uses arrays to store the offsets based on the topololgyId (with a special treatment for none).
The question is whether we want to merge this into the trunk, since the implementation is not well tested, yet. Additionally, about 5% performance in an explicit finite volume scheme might not be enough reason for some people. Opinions?
Attachments
Activity
Andreas Lauser (andreas.lauser@iws.uni-stuttgart.de) commented on the issue:
As a DUNE user, I think it would be great if this gets merged. This does not involve any user-visible API changes and it sounds like a vast performance improvement for the mapper. Did you any tests for the mapper itself without iterating over the grid? I would suppose that for a 5% improvement of the FV scheme, the speedup of the mapper itself must be much larger...
Robert recently was caught saying that the overhead of using dune compared to using ALUgrid directly is between 10 to 20% for an explicit finite volume scheme:
http://lists.dune-project.org/pipermail/dune/2010-November/007754.html
If your statement is correct, this patch potentially reduces it to 5 to 15%, which is pretty good IMHO :)
I also like the patch. The only thing I don't like is the 'magic' that computes the size of the offset vectors - but that's not you fault:
One reason to merge/unify topologyId and GeometryType was to be able to use vectors instead of maps as it's done in this patch. In order to allow this we should also export the maximal possible topologyId for each dimension. A static method int GeometryType::maxTopology(int dim) could do this. Perhaps one would also like an additional an enum in a template struct on GeometryType.
Similarly it would be nicer to use the index method of GeometryType (if we decide to add it) instead of using implementation details. Of course we'd like to know maxIndex() then.
BTW: Why not shift the indices by one or append an index for none instead of the more complicated pair<vector,int> approach?
I adapted the patch, see my personal git branch for this, and tested the patch. It looks good:
- I tested MCMGVertexLayout and MCMGElementLayout with an implicit one-phase Darcy flow from DuMuX using a cell-centered finite volume scheme and a box scheme with DoFs on vertices.
- The speedup is high, but I do not trust my data because I would expect less impact with our schemes. With a 300 x 300 YaspGrid the programs finish 10 % earlier (test_cc1p: 14.908s vs. 13.382s, test_box1p: 11.461s vs. 10.371s)
- I successfully ran the changes together with a custom MCMGLayout which uses codim=1 elements, but without measuring the performance.
What else must be done to include this improvement?
Wouldn't this be cleaner if it were based on the GlobalGeometryTypeIndex? I'd prefer having only a single implementation of the GeometryType -> index_offset mapping magic. Apart from that, GlobalGeometryTypeIndex would allow us to get rid of the second pointer lookup (we could replace std::vector<std::pair<std::vector,int> > with a plain std::vector without any additional work).
Finally (and slightly offtopic), shouldn't the Mapper reuse the index type of the underlying IndexSet? Hardcoding it to int seems kind of harsh…
-
I addressed Carsten's and Steffen's proposal and switched to GlobalGeometryTypeIndex. Please have a look before I'll merge my branch p/gruenich/fs854. I tested it again with cell-centered finite element method (element mapper) and box scheme (vertex mapper) from DuMuX. I compared the speed (dune-grid-howto/finitevolume.cc, int level = 9 * DGFGridInfo::refineStepsForHalf(), clang 3.3, dumux/optim.opts) and got another gain of approximately 5%: 1m23 current master 1m20 my branch with Martin's changes (96.4%) 1m15 my branch with my changes from today (90.4% / 93.8%) As always, I am skeptical concerning my benchmarking results.
As offtoppicly suggested by Steffen, I tried to use IndexSet::IndexType in the MCMGMapper but g++ didn't like it; yelled something about ambiguous method call. I didn't investigate it further (clang accepted the code), but I guess this change won't be that easy.
I had a quick look at your branch - looks good.
I think GCC might be upset about different method signatures in MCMGMapper and the base class (although I don't really understand why because there is no dynamic polymorphism). But still, in order to move to IndexSet::IndexType, I think you'd have to change the signatures in all four mapper classes at once, so that's probably a topic for a separate task.
One thing I was always wondering about regarding the MCMGMapper: The order of GeometryTypes within a given codim depends entirely on the grid (specifically the order of GeometryTypes in the vector returned by IndexSet::geomTypes() ). Is there some implicit rule that the GeometryTypes should be sorted according to their built-in default ordering, or is that unspecified? In PDELab, we actually moved away from that and use the global ordering implied by the GlobalGeometryTypeIndex. But that's mostly idle curiosity and shouldn't keep you from merging the patch ;-)
@Christoph: Go ahead. I did not commit the original patch, because we're not using the mappers in dune-fem. Being unable to test them, I did not want to break other people's code.
@Steffen: Your question addresses a lot of DUNE history.
Originally, not index was attached to GeometryTypes and their definition was rather random.
As some algorithms needed to find out what GeometryTypes might occur in a grid, the method geomTypes was added to the index set, returning a way of enumerating them. The use of std::vector is pure convenience (and should IMHO be replaced by an implementation-defined object returned by value). It does (to my knowledge) not imply any ordering.
The next step in history was the introduction of the generic reference elements. They are defined recursively, starting from the 0-dimensional object (a point) by building either a prism or a pyramid. Obviously, you can use bits to indicate which mechanism was used and, hence, you can assign a small number of indices to them. Apart from a small non-uniqueness issue and the type "none", this generates the GeometryTypeIndex.
In contrast to the geomTypes(), the ordering of GeometryTypeIndex cannot change during grid modification nor does is depend on the capabilities of the grid. However, being independent of the grid, it will cause some useless overhead for many simple grid implementations.
