Yes, you can encouter these problems while programming a computational code on unstructured grid. Year ago I worked on Navier-Stokes solver with immersed boundary method on unstructured grid. Any pointer to a node addressed a table with connections but I added many other topology infos loaded in the memory befor starting the run. A cache miss was still a possible problem but I found that a performance is generally increased at the cost of an increasing in the memory load.

Paolo Lampitella could provide many further infos

Hi Fillipo,

thanks for your prompt answer, this is what I had also in mind to have such a construct by increasing memory load ...

Aron

Hi Aron, actually I can't be of much help here. My experience is that introducing a Hilbert SFC ordering into a previously unordered code, without any particular optimization, can give a 30% speedup. Also, algorithms based on coupled set of equations (in contrast to segregated ones) can, under certain conditions, give some minor additional speedup. Basically, you vectorize on your local set of variables instead of the cells... but this might also turn out to be bad instead of good (a too large block of variables, like for gradients, might simply not fit the available cache in any case and if you need additional variables in the same loop... boom). I don't know what language you use for programming, but in Fortran don't use derived data types for this, just plain arrays.

Ciao Paolo,

mille grazie for this nice answer.

Yes, I can confirm that we get with Hilbert or other 10-30%. Another Yes, is about your comment on the coupled equations respectively systems of equation.

As for the programming language we r using Fortran but we are now thinking about porting it to C++ and using the latest features of this language, which give some advantages compared to Fortran.

Nevertheless, could u be so kind to explain why u would not use derived datatypes. I could e.g. by increasing memory usage and doing a construct like element(IE)%H(3) map everything into this derived datatype and acces than this stuff blockwise like element(1:BLOCK)%H(3) and work on the subset. As I understand your remark it seems that u anticipate some certain overhead by doing this and honestly I did this long time ago and it was not good but I am not sure that I really understand the reason for this. Maybe u are willing to comment also on that one.

Otherwise, Fillipo, Paolo many thanks u have been a great help since u underline my finding, which is not a solution to my problem but it shows that we r on the same line more or less.

Cheers

Aron

Ciao Aron,

lazy me for the incomplete sentence. What I meant is that, for example, for the 5 variables (p,u,v,w,T) of the coupled Navier-Stokes equations in 3D I would allocate an array like larger_type%var(5,ncells), instead of using something like larger_type(ncells)%var(5). So, in my opinion, use derived data types but only to contain large contiguous blocks (even if, in some cases, both approaches would produce the same code). Also, I might be wrong, but I'm not sure you can actually do something like element(1:BLOCK)%H(3). In any case, even if it is possible, the reason to avoid this is that you don't actually have control on the stride between element(ie)%H(3) and element(ie+1)%H(3). In some cases they would end up contiguous in memory but you're not ensured, especially if H itself is allocatable (I have experience of some cases with only one allocatable component in the derived type and still the compiler allocating an arbitrary amount of memory).

Obviously, all of this actually makes sense only if you effectively work with a coupled system (so that, for each face in a cell based FV framework, your flux function requires all 5 variables and their gradients from the neighbor cells). Otherwise this blocking only makes things worst.

P.S. Don't be so rude with Fortran ;)

The modern Fortran compilers totally rewrite the code when aggressive optimization options are used. We should better check (for example in the PGI or Intel documentation) if contiguous blocks are generated from the original source to prevent cache miss.

Hi Filippo,

I have checked all this using various profilers such as Vtune, valgrind etc. and there was significant chche misses and delayed instructions. What would be the benefit here of continues blocks?

As a side note on what Filippo wrote, in my block Gauss-Seidel iterative solver, I have higher performances if I write it down naively (i.e., accessing non changing variables within loops instead of pre-storing them upfront) and let the compiler do its job instead of trying to optimize myself. But this was kind of a corner case.

So, it might happen that you actually need to look at every specific algorithm and see if that specific piece of code is going good or bad and why. I mean, even looking at the specific line that causes cache misses and see if it can be cured independently from the data structure. In general, the less stuff is accessed in a loop the better it should perform (so that computing on the fly instead of storing sometimes becomes a valid, non obvious, option).

Also, the optimization level at which you try to optimize has a role. In my Roe solver, for example I have to multiply 4 5x5 matrices times a vector (why I do this instead of doing the product on paper a pasting it in the code is another story). When the product is written explicitly it goes much worst than matmul & co. at -O2, but slightly better at -O3.

Well, some experiments should be done to assess the relevance of the "ad hoc" coding when optimization options are used. For example, I used the PGI compilers when -fast option (associated to -O4 level) includes the -Mcache_aling option that aligns long objects on cache-line boundaries.

However, I worked on testing the performances a lot of time ago.

Hi All,

now after some time to come back. Solving the SWE we have 3 variables resulting of a subset of 3,3 mini matrices for each triangular element. As I have checked o the ordering ParMETIS is doing a good job in parallel and writing out the loop over the elements accessing the nodal values gives quite some dense pattern of memory access, which is nice. Still, the overhead of the "data extraction process" of the elemental values of the state variables at the nodes is a killer.

So, the next aim is to look on SIMD stuff and try to code everything "elemental" mini functions and overload them to increase vectorization ... let's see.

Cheers

Aron