literature on scientific computing

Question

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.

Answer 1

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

Answer 2

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