Developing finite volume (FVM) code in C . General advice

Question

I need to develop a FVM code in C (The multiscale FVM method for heterogeneous media). I know that:

• Only uniform rectangular grids will be considered (2d now, later 3d)

• Sparse systems will be large

• Need to parallelize my code eventually (will be solving local problems on fine grid and use them to solve PDE on coarse grid)

• Speed is key

I used Matlab exclusively during my PhD. I am now faced with the challenge of C programming which is new to me.

I fully intend to use existing libraries as much as possible, but I am quite puzzled about my options and to what extent I should use them.

1. Must I write my own program and assemble the matrices and then look for solver libraries, or are there alternatives ?
2. Which library would you recommend for this? Is PETSC my best choice ? (I couldnt find many sparse solvers in C)
3. I looked at GSL, LApack, etc and my question is does using a library entail using only the data types/headers that they offer ? For instance GSL has a vector and matrix header file which seems to have its own data type.
4. Even simple steps, like setting a matrix A= B requires several lines of code in C, is there any library that offers some of these basic "matlab" luxuries, basically readymade functions in C for such tasks ?

I realize some of my questions might be very elementary, considering my lack of familiarity with C and the non-adhoc way to program numerical methods. I apologize.

I understand the question is subjective in nature, but your answers will point me in the right direction.

Answer 1

Some thoughts from someone who has worked a fair amount in compiled languages, and has done a tiny bit of FVM:

1. Typically, if you have experience programming in C, you sketch out a high-level description (pseudocode) of what you would like to do. Then you look for libraries that might implement the data structures and capabilities you need for your high-level description. Then you learn how to use those libraries by writing simple programs that do simpler versions of what you want to do. Sometimes, libraries have example programs. Use them. If you can't find a library that does what you want, you write code to do it yourself. Then, once you have some handle on what you want to do, and how you use the libraries you want to use, you write your program. This process may be iterative; for instance, you may find that you neglected some crucial step in your high-level description, and then you have to find a library that implements the capabilities you need for that step, and repeat the process. Write programs to solve small problems first, then work up to bigger problems. Trying to write the program that will solve the big problem in the first attempt is just asking for problems.

2. From your description, PETSc is probably a reasonable choice. It's well-documented, the developers are responsive and helpful, and it has many examples. PETSc's certainly not your only choice, but it's probably the most convenient library that implements standard numerical methods from time steppers to linear solvers, and it makes it relatively straightforward to write a code to run in parallel (assuming you have some basic idea of how to do so).

3. If you're asking this question, you should read a book on C and do the exercises. The book by Kernighan and Ritchie is good, as is Learn C The Hard Way, by Zed Shaw. Also, read a book like Code Complete by Steve McConnell, or Clean Code by Robert Martin, so you learn decent programming style, and so you learn to write something that other people can read.

   Typically, you include the library header files and link any static or dynamic library files accompanying those headers. Libraries provide an Application Programming Interface (API), which includes data types and functions (typically operating on those data types). As you program more, you'll see what I mean.

4. Not really, no. In using C, you're trading off programmer productivity and the convenience of automatic memory management for execution speed and manual memory management. If you want something that's more MATLAB-like, use Python. There is a Python interface for PETSc called petsc4py. It is not as well-documented as the PETSc library. For finite volume methods, if you're willing to live with explicit time integrators, you could try PyClaw, also written in Python. Pure Python programs will be slower than pure C programs; however, many Python libraries write performance-critical portions of their code in C, so you get most of the performance benefits of C along with the convenience and productivity benefits of using Python. (You could do something similar in MATLAB with MEX files.)

Answer 2

DEVSIM is a finite volume code that I've implemented in C++ with a Python interface. It uses a generalized PDE solver, meaning you can easily implement the Poisson equation in it (Please see the capacitor example in the program documentation). It can create a rectangular grid (as triangles) or read in a mesh using created in gmsh. You are welcome to try it out to see if it suits your needs.

If you are restricted to C (not C++), you could try UMFPACK or SuperLU solvers for doing the sparse matrix calculations with your own implementation. If you can use C++, you can use PETSC, which could use these libraries as underlying solvers.