# Effecient CFD programming techniques

I'm trying to make highly efficient CFD programming complex for solving combustion problems. I've finished writing core which realises mathematical model, and now I'm concerned about code performance. Combustion problems are very computationally expensive, with no doubt I'll have to make parallelization of my code, but for now I want to make highly efficient serial code, to use it as a base for parallel program.

For instance there is simplified version of what is going on in program:

I have the main program where time loop is located.

PROGRAM Flame
USE Geometry
CALL ALLOCATION  <- Procedure there memory is allocated for X array
DO while ...
Time=Time + TH
CALL CALCULATIONS
END DO
END PROGRAM


Modules with all global variables declarations

MODULE Geometry
double precision , allocatable :: X(:) <-- one dimensional array for computational grid.
integer, dimension(NTR) :: KDIF <-- has a mean of pointers on some variable in cell (pressure for example)
END MODULE


Subroutines for calculations and data processing.

SUBROUTINE CALCULATIONS
USE Geometry
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
IJK = 0
DO I=1,NX <- number of cells in one dimension (such cycles for all dimensions)
DO I1=1,NTR
X(IJK+KDIF(I1)) = ... some floating-point operations <-- Possible hot spot
I1 = I1+1
END DO
X(IJK+8)=X(IJK+3)-A/X(IJ+5)-(X(IJ+10)-B)*C/(D*X(IJ+5)) (for example)
IJK = IJK + NKX  <-- place between two cells in
END DO
END SUBROUTINE


There is much more subroutines and modules in program, but they all have structure similar to that I gave as an example. I've tried to profile my program with VTune Amplifier, and located some hot spots in places like I pointed in code. In Assembly of hot spots like this most time consuming operations are generally movesd qword; move ebx, dword; etc. I'm not familiar with assembly codes but I think it is related to memory operations. Actually I have no idea how to fix that hot spots.

So there are some questions that I need to clarify. Is that style of coding modern? I mean maybe I should use another techniques for such programs, like pointers, FORALL and WHERE operations, which can be faster and more memory efficient, more appropriate for CFD problems? Is such program structure efficient (modules for data, subroutines, allocatable arrays, one-dimensional array even for multidimansional problems)?

• I would recommend looking into open source packages to get an idea of how they are put together for performance. OpenFOAM is an open source CFD solver (I have personally never used it, so I don't know what combustion/chemistry it would support) that would probably be a good place to start. Optimizing at the assembly level is likely a bit premature until you ensure you have chosen the correct solvers and algorithms for your problem. – Godric Seer Apr 13 '13 at 23:57
• I have had experience with OpenFOAM package. It has quite complex computational core written in C++, although I think it is good idea to try to look there. Algorithm is fully fit to my tasks (Euler-Lagrange Large particle method), and was tested for combustion problems, so I do not think about changing the algorithm. My task for now is to create its highly efficient implementation in Fortran. – Yakovenko Ivan Apr 14 '13 at 9:00
• It seems that you are using explicit time-stepping. Maybe you can design an updating ordering of dof's to have cache-efficiency? – Hui Zhang Apr 14 '13 at 10:54
• Could you please give a simple example of updating reordering? I was thinking about structure of my main array X, whether such implementation of computational grid is effective. May be I should reorder it so data needed for calculations will be closer in array, now it is ordered cell wise (IJK(cell index) + NKZ(Z neighbor)+NKY(Y neighbor)+NKX (X neighbor)+V (Variable index)), may be it is more efficient to order it variable wise (V+1(Variable V in Z neighbor)+NZ(Variable V in Y neighbor)+NZNY(Variable V in X neighbor)*. But I don't quite understand how to make updating reordering... – Yakovenko Ivan Apr 14 '13 at 13:15
• See e.g. the reverse cuthill-mckee en.wikipedia.org/wiki/Cuthill%E2%80%93McKee_algorithm I believe there is a Fortran implementation. – Reid.Atcheson Apr 14 '13 at 21:55

Optimizing code is a very broad problem. Before you get too far into the details of an implementation, you should make sure your models and algorithms are doing what you want them to as efficiently as possible. A good algorithm is better than a good implementation and premature optimization is the root of all evil.

That much being said, here are some comments on the code you've given:

• Don't use implicit typing, ever; it turns typos into bugs.
• There are mixed opinions on forall, where, etc. This discussion is particularly interesting. Summary: forall, where, etc don't provide much, if any, speed, but they can be nice for representing algorithms.
• Don't pass array slices as arguments; they can cause unintended temporary allocations and copies.
• allocatables are fine; only switch to pointers if you need more functionality
• Modules are fine. I personally like to reserve modules for parameters or variables that are nearly parameters (and routines), and move data around explicitly. Fortran is pass by reference, so large array arguments are no more expensive than scalers. I've seen successful codes do it either way.
• Thanks alot for your answer, it was very useful! I think that after chosing numerical algorithm (Euler-Lagrange Large particle method) for my problem subsequent stages are clearly defined. In according with method algorithm I should translate mathematical formulas into programming language. That is also done, but I think that it may be performed better in terms of optimizing. So now I'm searching for any programming techniques that can make my code better to be sure, that i have no optimization issues on serial level before i begin code parallelization. – Yakovenko Ivan Apr 14 '13 at 12:45
• Keep in mind that parallelization will likely require some changes to the data layout and flow of the program. Those changes can move your hot-spots' around and dictate which optimizations are efficient. Some of the effort you spend now on the serial code will likely end up going to waste... – Max Hutchinson Apr 14 '13 at 15:58
• forall is rarely helpful for performance because it has semantics of global vectorization when evaluating the right hand side, which often requires temporaries and is pessimal for cache reuse (far more important most of the time; you can do 50 flops in the time required for one perfectly-prefetched streaming load). If the compiler can determine that the RHS evaluation and assignment is safe, it will generate the same code in both cases. In general, it is easier to optimize nested do loops. – Jed Brown Apr 14 '13 at 16:06
• Array slices are often bad for performance as well, especially if you ever take a non-contiguous slice. Calling a function with such a slice involves allocating and filling a temporary, calling the function, then unpacking the result (whether or not it was modified because the caller doesn't know). – Jed Brown Apr 14 '13 at 16:08
• @MaxHutchinson The first search results all agree with my own experiments from about 10 years ago that forall is slow. It wasn't really intended for this, and its semantics are not made for cache locality. Some people have the misconception that forall comes with the equivalent of IBM's independent` extension, at least when reasoning about performance. compgroups.net/comp.lang.fortran/forall-vs-do.enddo-vs-where/… stackoverflow.com/questions/4122099 – Jed Brown Apr 14 '13 at 17:46