Object-oriented programming on finite difference method

Question

generally, it is natural to use procedural programming approach (PP) to solve a partial differential equation by finite difference method (FDM). That is, one (1) defines matrices to store the properties at nodes and fluxes in between the nodes, (2) assembles linear systems based on discretization to the controlling equation , (3) solves linear equations and (4) moves into next time step. As each controlling equation solves one property across the whole domain, defining a matrix of one property is very convenient.

but so far I realized some (maybe) disadvantages of using this methods, which is that the modules (or subroutines) used for process (2) do not easy to be ported to another program solving another fundamental equation. basically if one wants to solving another problem, it is required to start from scratch (except for matrix solver of course).

I have been adviced to use object-oriented programming method (OOP) to handle the problem, which is very convenient to modulize the code, and the modules are quite reusable.

for my understanding, by OOP method, one may not define a matrix storing the property for the whole domain, but define each node with many properties (using objects). This method is quite robust in some Lagrangian methods like smooth particle hydrodynamics and lattice Boltzmann method. but this method may not quite convenient to solve linear system as the matrix has been disassembled.

In contrast, Nowadays many programming starts to support OOP, even fortran starts to support objects. In addition, currently there are many FDM, Finite element method (FEM), Finite volume method (FVM) packages using OOP concept (e.g. openfoam), one only requires to define the fundamental equations and then desensitization seems to be done using the same codes. I think that may not be easy to use PP to handle.I must have something missing about OOP ideas.

of course this problem may not limited to use finite different method, but also to use FVM and FEM which to me, seems more difficult to construct by OOP concept.

any ideas? thanks in advance!

Ps: generlly I use fortran, matlab, R, python deal with CFD and multi-phase flow problems in porous media using PP.

Answer 1

OOP has great utility in scientific programming, but the natural way one would write an "OOP" scientific program after education in non-scientific software engineering may be a bad way to tackle the problem. Anecdote:

• My first semi-serious CFD code as an undergrad, I wrote a finite-volume 2D Euler equation solver with excessive use of OOP. Intuitively, I figured a "cell" would be a good class, each with pointers to "face" classes that are shared amongst cells. Then I iterated through cell and face objects. i.e., I had a few million tiny objects of "cell" and "face" type. It was incredibly slow -- every access of these elements resulted in tons of dereferenced pointers.

As a general rule of thumb, I'd suggest taking the advice of the CUDA people (and for the same reasons -- it's not news to people who know about high performance code):

• Small structures of large arrays are better than Large arrays of small structures.

Check out Damian Rouson's "Scientific Software Design: The Object Oriented Way" for a good take on ways to use OOP for high performance code. I have some serious issues with details in there, but it's still pretty good.

Answer 2

Absolutely nothing is "natural" about using procedural programming for finite differences or finite elements. In fact, all widely used software packages for finite element methods that have come up in the past 15 years are object-oriented (e.g., libmesh, openfoam, my own library deal.II, fenics, ...) and the same is true for linear algebra packages (e.g., petsc, trilinos). That doesn't mean that every single thing (like a cell) should be an object, but certainly that things like meshes, matrices, finite elements, are all implemented as classes. There is absolutely no question in our community that this is the only possible way to deal with large-scale software.

There are several other questions on the SciComp StackEchange site that deal with similar issues.