3
$\begingroup$

Maybe, this is a dump question. But anyway.

I developed a finite difference code for a class of problems in dynamo theory. I used GNU Octave (MATLAB) which is good for testing. The problem size scales quadratic. Hence, the number of unknown becomes relatively large.

I am considering to rewrite my program in a different programming language. In my impression the most common languages used in scientific computing are FORTRAN, C++ or Python. The point I would like to make is that I don't want to reinvent the wheel.

Some aspects to get the idea:

  • I would like to use existing libraries/ packages such as LAPACK, ARPACK, UMFPACK, PETSc,... for linear algebra applications. I don't want to write Gaussian elimination again.
  • I don't want to program every matrix multiply in for-loop.
  • I would like to export unstructured data to the vtk file format. I would like to use 2-D plotting tool like in GNU Octave.
  • I would like to use MATLAB like methods such as sparse or blkdiag.
  • The possibility of running on multiple processors would be appreciated, i.e. openMPI.

I am neither looking for a book such as Numerical Recipes, a book on numerics of PDEs nor an introductory textbooks on programming. I am looking for a book that describes the general approach. How to structure the program? How to make us of existing packages? An example, such as the standard 2-D Poisson equation would be good.

|cite|improve this question|||||
$\endgroup$
  • $\begingroup$ You can try PETSc tutorials. $\endgroup$ – stali Feb 24 '15 at 13:44
  • 1
    $\begingroup$ How many equations are there in the model you would ultimately like to solve? I am doubtful the MATLAB (or Octave) sparse matrix solvers are to blame for the poor performance. Have you used the MATLAB profiler to determine where the time is being spent? Even if you end up rewriting your code in another language, this exercise will likely be useful. And I wouldn't automatically assume it is impossible to create a MATLAB implementation with acceptable performance. Sometimes you can substantially improve the performance of a MATLAB program by rewriting a small portion as a MEX file. $\endgroup$ – Bill Greene Feb 24 '15 at 14:13
  • 1
    $\begingroup$ @stali, and/or the deal.II and libmesh tutorials. $\endgroup$ – Bill Barth Feb 24 '15 at 14:32
  • 2
    $\begingroup$ It's not quite clear what exactly you are looking for here -- your question mixes language requirements, specific problem details and a book recommendation. Can you edit your question to be more explicit about the kind of answer you are looking for? (It's OK to give some background, but the actual question should be front and center.) $\endgroup$ – Christian Clason Feb 25 '15 at 14:45
  • 2
    $\begingroup$ You should also check to see if the answers to this question help you: scicomp.stackexchange.com/q/10563. I think the literature there is exactly what you're looking for. (I'm not aware of any books apart from Hans Petter Langtangen's, which are mostly Python-specific.) $\endgroup$ – Christian Clason Feb 25 '15 at 14:52
5
$\begingroup$

Honestly, the best ideas have already been said. I'll try to synthesize my thoughts anyway.

First, the best way to write a program is whatever way gives you the results you need, for the least development time. For some applications, you NEED bare-metal performance, but a lot of the time you don't. Since you've been working in Octave up until now, you may not fall into the bare-metal category, unless you're suddenly increasing your problem size dramatically. Even if you are, the C vs. Fortran speed debates are mostly immaterial at this scale. Even a factor of 10 doesn't matter much if it's 0.1 second vs 1 second.

Second, my experience has been that most of what a code does is best handled by a high-level language such as Matlab, Python, etc. There are usually only a few small segments of code that really benefit from the extra programming effort to really optimize them. Often, these are things like inner-loop solvers, complicated compiled functions, and so on. It can be very beneficial to write the majority of a program in a high-level language, and then call low-level library functions as-needed for the really performance-intensive stuff. .mex files are Matlab's way of handling this (it looks like Octave has a clone of this, too). Python is especially good at it, too, with things like f2py, Cython, numba, and so on.

Third, and finally, the most important aspect of this by far is the software ecosystem that will be available to you. Matlab isn't widely used because it's so blazing fast; it's widely used because it can plot data, perform LU-decomposition, numerically integrate, and much more, all with simple commands, and knocking a few months (or years!) off of development time is almost always worth the performance hit.

To sum up, you need to quantify your needs (what's your problem size, and what's the acceptable length of time for a solution), examine what you already have (for re-use and/or porting), and make a list of your ecosystem requirements (vtk, linear algebra, etc.). Depending on the specific answers you come up with for these, there are many possible "best" options. Some of these might include:

  • Keep your Octave code, but refactor it to use external libraries using .mex files where the performance requirements justify it.

  • Rewrite your code in Python, for access to a wide variety of libraries and excellent mixed-language programming tools. This is what I have done, fairly successfully, in my work.

  • Write a new, completely parallelized HPC code using PETSc, MPI, etc. Unless the problem size dictates that you run this code on a supercomputer, my opinion is that you should stay away from this one. A friend of mine migrated a serial CFD code to PETSc once, and it took him years. YMMV

|cite|improve this answer|||||
$\endgroup$
3
$\begingroup$

I think you are going to have to choose what your primary criteria is: performance, modularity, or portability.

I took a highly optimized Fortran 90 code and tried to modularize it and make it more object-oriented, and the performance became 3X worse in the object-oriented version. I found the best I could do to have the best of both worlds was to make OO data structures that were all pointers pointing to the original non-OO data structures. Fortran uses highly efficient array strides which makes it very fast for doing linear algebra computations. The problem is that Object-oriented (OO) approaches destroy the optimal cache coherence you get from loops in non-OO approaches. There was a pretty good discussion on that here.

If you are primarily concerned about performance, I would use a combined Python / Fortran 90 approach and connect them via the f2py module. Fortran 90 in general is faster than C or C++ especially for linear algebra (although depends on the platform), and has good support for many of the numerical libs such as OpenMPI, Lapack, etc. I would do all the compute intensive (e.g. linear algebra) stuff in Fortran 90, and do all the coding intensive stuff in Python (e.g. vtk). There is an example of solving Laplace's equation $\nabla^2u=0$ using a combined Fortran 90 / Python approach via f2py here. Python has one of the largest scientific computing communities for a modern high-level language, and has good support for a number of modules, especially numpy, scipy, and also the vtk module. There is an example of reading unstructured grid via Python's vtk module here.

On the other hand, if you want to be more modular and/or portable, you may want to look into using Java, but you will have to sacrifice performance. I haven't been too impressed with Java-based simulations especially in terms of performance.

Concerning reference material, I don't think there are any books covering the software architecture of these kinds of problems. However, you may want to checkout this paper titled "Using design patterns in object-oriented finite element programming".

I don't have any experience with it, but there is a rapidly growing interest in coding these kinds of problems in Julia. So, you might check out this as well: http://www.codeproject.com/Articles/579983/Finite-Element-programming-in-Julia

|cite|improve this answer|||||
$\endgroup$
  • 1
    $\begingroup$ Your criticism of OOP, isn't that normally known as structure-of-arrays vs array-of-structures? That's not really a problem with object-orientedness, strictly speaking, because it doesn't use any OO features, like polymorphism or encapsulation. Instead, it's purely a memory layout issue. $\endgroup$ – Kirill Feb 26 '15 at 16:30
  • 1
    $\begingroup$ Regarding your commentary on F90 vs C/C++: I think there's a lot of computational science lore that gets spread by word of mouth that F90 is faster than C/C++. In practice, that's not really borne out; Fortran 77 didn't have to worry about pointer aliasing, but these days, if you know what you're doing, all three languages are competitive. If anything, I'd argue that Java is less portable, since its MPI bindings are crap. $\endgroup$ – Geoff Oxberry Feb 27 '15 at 0:32
  • $\begingroup$ Python is mostly used as syntactic sugar and I would stay away from it. A recent benchmark (advective transport) of Python, C++ and Fortran is at downloads.hindawi.com/journals/sp/2014/870146.pdf $\endgroup$ – stali Feb 27 '15 at 3:09
  • $\begingroup$ Another benchmark: julialang.org/benchmarks $\endgroup$ – Wes Feb 27 '15 at 5:15
  • $\begingroup$ @Wes Good insights about performance aspects. I guess I will have to do more research about the languages and libraries. So far, I tend to use python, i.e. numpy, or the PETSc library. $\endgroup$ – sebastian_g Feb 28 '15 at 9:12

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.