Data structures for finite volume code: Arrays vs Classes

Question

I have to write a finite volume code for Magnetohydrodynamics (MHD). I have written numerical code before but not at this scale. I just wanted to ask which will be a good choice, using a data structure (object orientated approach) with classes or just using multiple arrays for different properties, in terms of speed, scalability etc. I plan to write the code in python, and use fortran for numerically intensive part.

An example for class in python would be

class Cell:
   def __init__(self, x, y, z, U):

Arrays can be simply defined as

x[nx][ny][nz]
y[nx][ny][nz]
z[nx][ny][nz]
U[nx][ny][nz]

etc.

Answer 1

Simple answer: in modern python every data type is a class, so formally there is no difference between the two solutions you proposed. (Please remember to use new-style classes: classic classes are obsolete! See http://docs.python.org/2/reference/datamodel.html#new-style-and-classic-classes)

Now the question should be: how do I organize an efficient data structure in python? There is no doubt that the very idea of organizing the cells as an array of class Cell instances is way too inefficient. You will end up with a mess of pointers and non-contigous data organized like a complicated linked list. You have of course the ability of easily insert new cells in your list: but do you need this feature? On the contrary you will have non contiguous data storage, and you have to access every cell by different levels of indirection.

If you organize your data as a numpy.ndarray then data is memory-contiguous, and accessing different cells is simply done striding through your memory block: space efficient (no memory wasted for pointers) and fast.

As pointed out by Ethan, OO concepts should be used, but at higher level, once an efficient low level data structure has been implemented, usually through numpy.ndarray's.

OO programming means binding data to the methods that operate on the data itself at higher level of abstraction. (An example: I implemented a FEM code in which the stiffness matrix was defined as a class with a method for sparse super-nodal cholesky factorization. The first implementation was in-core: when an out-of-core implementation was needed, this was obtained via inheritance and minimal adjustments to the underlining data storage. Almost 100% of the super-nodal cholesky code was reused.)

A last comment, but crucial: an efficient numerical procedure is the result of a smart mapping of an algorithm and a data structure to your target computing architecture. If you start with the wrong data structure, there is no way of recovering efficiency, without a complete rewrite.

Answer 2

I was pondering this a few days ago (also in Python). Personally I don't think that object oriented programming is always a good fit for numerical programming. You can get distracted with designing the classes rather than just solving the equations. I prefer to stay with simple functions, and with numpy you can have your equations vectorised so the number of lines you need is very few. Numpy is pretty fast because the actual computations are done with a C (or FORTRAN?) back end.

What I would recommend you do,

1. Write a Python script that solves the simplest possible version of your problem using a functional approach with numpy. For example, have everything in arbitrary unit and try 1D (or 2D) only. It is perfectly fine at this stage if the code is messy. The important thing is that you are moving forward with your project.
2. Once you have something that works. Identify where the code is verbose and refractor. At this stage you can play around with different ideas on how to simplify your code. Maybe introduce functions where you notice you are repeating yourself. You can compare with the original version so you know that you are not introducing bugs.
3. Decide if object orientated approach will reduce the complexity of the code further.

The main message is don't start writing classes until you have already solved the problem in the simplest possible way. It is only through gaining the experience of solving a problem you will know how to define your object orientated interface. If you do this before hand it is likely to just get in the way.