Comparing various implementations/software packages for large-scale finite element simulations

Question

I currently use FEniCS and Deal.II to solve various FEM problems. I am also writing my own implementation of these problems by directly implementing the data structures, routines, and solvers within PETSc. What kind of comparisons can I draw between these three various implementations given the same FEM problem/discretization? Particularly in the context of high performance and parallel computing.

Naturally, the first thing to do is check the numerical accuracy - do I get the correct solutions. Then I could check the strong scaling, weak scaling, and associated wall-clock time across multiple processing cores. But what else can I do? For instance, is there a way to determine how "efficient" these frameworks are algorithmically and in terms of machine hardware?

My goal is to show people which of these implementations will yield the fastest possible time given a specific problem size and FE discretization, but I think it would also be interesting to show how "credible" these timings are. For example:

1) Why does implementation X take 45 seconds to solve a problem but implementations Y and Z only take 30 and 25 seconds respectively?

2) Why can't any of these three take only 4.5 seconds?

3) How much of the time within the computational framework is spent doing useful work versus waiting on/accessing memory/cache?

Sure memory bandwidth is likely the answer to why these implementations don't achieve the ideal or perfect time, and benchmarks like STREAM can tell you what you can expect given your machine. I am wondering if anyone has a methodology, tool, or suggestion to properly determine or quantify these.

Thanks,

Answer 1

Unfortunately there's no tool for this. You can run each on a variety of input sizes to establish the computational complexity they appear to have, i.e. the $f$ in the $O(f(n))$ that characterizes each code. This can point you towards what underlying algorithms each is using and verify or not that something that should be $O(n)$ is actually implemented to achieve this. You can run fixed $n$ and vary the number of processors $p$ in order to try to determine the efficiency of their parallel implementations. And then you run a fixed problem on a variety of CPU frequencies and memory frequencies to try to determine how sensitive they are to hardware characteristics.

These studies can help you lay out a handful of mostly orthogonal axes in the design space. This space is actually much bigger and quite vast, and it doesn't even begin to take in the various input parameters to the solvers (max iterations, stopping criteria, algorithm selection, etc).

As best I can tell, basically no one does this. You usually figure out pretty quickly whether you have accidentally implemented something $O(n^2)$ when you meant to be $O(n)$. With practice, you can look at the code to most core algorithms and be sure you got the complexity right, so in the end, when comparing two implementations that are supposed to have the same computational complexity, you're trying to compute whether they have achieved through implementation the multiplicative and additive constants the theory predicts. This can be very tedious.

Once you've done all of this, then you have to spend the time comparing the implementation, the compilation to machine code, and the hardware itself to see why you didn't get anywhere near peak performance. People spend years doing this.

Answer 2

As one of the library's authors, I would of course love for deal.II to come out on top with this comparison. But I suspect it may not, and the answer lies in a factor you omit from your comparison: how long it actually takes to implement your code.

Few people in academia with the skills to implement a FEM code from scratch spend more time solving PDEs than actually implementing them. In other words, the time to implement things is more valuable then the actual run time of your code. This is what brought us libraries such as deal.II, FEniCS, PETSc, etc: Yes, they may not be the very most efficient implementation of whatever it is you want to do, but they are generic in the sense that they can be used for a large variety of problems, and they are comprehensive in that they get you 90% there without having to write very much code. For example, deal.II's parallel Laplace solver, step-40, has only around 135 lines of code (semicolons). That's likely 2 orders of magnitude less than you had to implement if you wanted to do it from scratch -- saving you from writing 20,000 lines of code yourself (about a year's worth of work). Furthermore, switching to a higher order element is a one line change. FEniCS may even get you down to 20 lines of code, though the saved 115 lines of course don't quite add up to very much compared to the initial savings.

Point being, performance and run time are only one metric. They matter if you intend to burn hundreds of thousands or more CPU hours. For the rest of us, development time is more valuable.