Optimize data structures for performance without losing readability/flexibility.

[pcarruscag]

Having spent some time on implementation details of linear algebra topics, and a bit of time on small algorithmic improvements to a few numerics classes, I would now like to focus on another performance bottleneck we have. One that always seemed quite daunting to tackle (to me) due to the amount of changes it would require, but recently I had a silly idea to reduce the invasiveness of the "procedure".

Discontiguous data structures are absolute performance killers for low order FVM.
(skip this paragraph if this is obvious)
What data structures? All the small ones the solver needs to "visit" to compute residuals.
Why? Mostly because they destroy locality. The useful data stored by consecutive nodes is very far apart. Modern CPU's try to hide the latency of accessing main memory (RAM) by caching frequently used areas of that memory, this effort is very easily undermined by accessing memory at random or predictably but at very large strides.
Let me give you some numbers, for CNSVariable in 2D the average distance between the primitives at consecutive nodes is 200 doubles or 1600 bytes. "What?! is the compiler dumb?" No, virtual methods take up quite a bit of space sizeof(CNSVariable) is 424 bytes, before it allocates any data.
And that can only get worse as more functionality is added to the code.

This is one of those "well known" things at least in a qualitative way, but let us try and quantify the cost using an approximation of our typical use case scenario.
Spoiler alert: It is up to one order of magnitude, if this does not surprise you feel free to go straight to "Proposed solution".

We will have a variable class that stores some data with some overhead (the 200 doubles):

class VariableBase
{
  protected:
    double* primary;
  public:
    VariableBase(int nVar, int overhead);
    virtual ~VariableBase();
    inline double* GetPrimary() {return primary;}
    inline virtual double* GetSecondary() {return nullptr;}
};

VariableBase::VariableBase(int nVar, int overhead)
{
    primary = new double [nVar+overhead];
    for(int i=0; i<nVar; ++i) primary[i] = 1.0;
}

All this does is give access to some primary variables held by the base class and allow derived classes to store their own secondary variables, e.g.:

class VariableDerived : public VariableBase
{
  private:
    double* secondary;
  public:
    VariableDerived(int nVar, int overhead);
    ~VariableDerived();
    inline double* GetSecondary() {return secondary;}
};

Then lets say we have dummy numerics classes.

class NumericsBase
{
  public:
    virtual ~NumericsBase();
    inline virtual double ComputeResidual(const double* a, const double* b) {return 0.0;}
};

class NumericsDerived : public NumericsBase
{
  private:
    int n;
  public:
    NumericsDerived(int nVar);
    double ComputeResidual(const double* a, const double* b);
};

NumericsDerived simply computes the dot product of a and b.
The use case is to loop over edges that connect nodes a and b passing the data of those nodes to ComputeResidual and using the result to update an array.

VariableBase** node = new VariableBase* [nNode];
for(int i=0; i<nNode; ++i)
    node[i] = new VariableDerived(nVar,overhead);
...
for(int i=0; i<nEdge; ++i)
{
    int p0 = connectivity[i].first;
    int p1 = connectivity[i].second;

    double r = numerics->ComputeResidual(node[p0]->GetPrimary(), node[p1]->GetPrimary());

    residual[p0] += r;
    residual[p1] -= r;
}

I hope that looks familiar.
connectivity is a vector<pair<int,int> > so nicely contiguous (in SU2 this information comes from an array of CEdge*). Note that we are never talking about sequential access here, this level of indirection will always be present.
To create a somewhat realistic access pattern we create 2 edges per node connecting it to nodes +- some distance away, this distance is small for end nodes and largest for the middle node (I assumed the maximum distance to be sqrt(nNode)). This should vaguely resemble an RCM ordered 2D connectivity.
We will measure the cost of this edge loop relative to using a row-major storage of our variables.

double* solution = new double [nNode*nVar];
for(int i=0; i<nNode*nVar; ++i) solution[i] = 1.0;
...
    double r = numerics->ComputeResidual(&solution[p0*nVar], &solution[p1*nVar]);
...

OK, so we run this for 1000000 nodes with 4 variables per node and a overhead (or "space between variables") of 200 and we see that CVariable is 8.4 times slower.
But we can actually confuse the CPU further by accessing the variables through the virtual method, i.e.

...
    double r = numerics->ComputeResidual(node[p0]->GetSecondary(), node[p1]->GetSecondary());
...

and now we measure our loop to be 10 times slower.
If we remove all overhead the loop is only 1.5 and 1.9 times slower respectively, it is really the enormous distance that kills performance.

Proposed solution
Now, I would like to get better performance, but I do not want to re-write the entire code, and I recognize the usefulness of polymorphism and that encapsulating data keeps the code more manageable and readable.
To that end I propose that CVariable should allocate data for all nodes, so CSolver would have one CVariable and not one per node e.g.

class DenseVariableBase
{
  protected:
    double* primary;
    int n, i;
  public:
    DenseVariableBase(int nVar, int nNode)
    {
        primary = new double [nVar*nNode];
        ...
    }
    virtual ~DenseVariableBase() {delete [] primary;}
    inline DenseVariableBase* operator[] (int iNode)
    {
        i = iNode;
        return this;
    }
    inline double* GetPrimary() {return &primary[i*n];}
    inline virtual double* GetSecondary() {return nullptr;}
};

To keep the same semantics of accessing variables node[iNode]->Get() (and thus avoid changing the entire code) the class now keeps track of what node we are accessing (member variable i) and all its methods have this in consideration see for example GetPrimary above.
We tell the class what node we want via operator [] which returns a pointer to the object itself so we can call other methods on top of [].
This does not exactly recover the same interface and it can create some bugs because one would be able to call methods without knowing what index they are on (it is also not thread safe but that is not a problem now).
So the final piece to recover the interface and shield direct calls to methods without specifying the index is a lightweight wrapper.

class DenseWrapper
{
  private:
    DenseVariableBase* nodes;
  public:
    DenseWrapper(int nVar, int nNodes);
    ~DenseWrapper();
    inline DenseVariableBase* operator[] (int iNode)
    {
        return (*nodes)[iNode];
    }
};

With this there is no way of getting to the variables without going through [] first.
If we run this code we see that it is only 1.08 and 1.57 times slower for the non-virtual and virtual accessors respectively.
This could be an intermediate solution while the CVariable methods would be migrated to allow passing the index as argument e.g. field->Get(iNode) and the code is updated gradually.

However, I think the wrapper could provide a way of dealing with the loss of versatility that contiguous storage creates...
There are no free lunches and this performance improvement would come at the cost of losing the ability to have different types of variables in different parts of the domain. I think currently only the elasticity solver uses this to keep "boundary" and "internal" variables.
If this versatility is really required the wrapper can keep special variables separately, that it would know to be special via a boolean vector and know how to access via a hash table (for example).

Alternatives
We could make the storage continuous by using a memory pool i.e. allocate outside the class and give the class a pointer to its memory location. This is still not truly contiguous as the pointers the class holds are variables themselves and so they would still be separated (albeit by a smaller huge distance), thus the amount of overhead would still influence performance.

Right, that is enough writing... please let me know what you think.

[talbring]

thanks @pcarruscag for this explanation and suggestions. I think most of us are aware that the virtual functions cause an overhead in runtime. Still, never thought it would be that drastic. However, I would be interested to see on how much faster it is in practice, since most of the work is currently done in the Compute_Residual of the numerics classes and the solution of the linear system.

But I like your proposed solution as it really keeps the code readable. So I don't really see a problem in adopting that in the code.

[pcarruscag]

Yes, I did not make that clear, the linear solver fraction of the time cannot be accelerated by this. But everything else can, from gradient/limiter computations to the Compute_Residual functions, as all those need to wait for data, it is however not very easy to measure how long that wait is compared to the rest of the computations.
In the 2015 joint work between Stanford and Intel they reported a 1.5x speed-up from this type of change for a case where the linear solver used 24% of the time. I do not know how heavy the CVariable infrastructure is now compared to then... We will see :)

[rsanfer]

Hi @pcarruscag,
thank you very much for the insight and the well thought proposed solution. I wasn't aware that there was such a loss in performance, it seems clear that something needs to be done.

The only question I have is regarding the following comment

There are no free lunches and this performance improvement would come at the cost of losing the ability to have different types of variables in different parts of the domain.

I'm using this feature more and more, I was planning to extend it in the structural solver to remove the if fsi statements, and I also used it in the rework of the mesh solver, as defining variables at the interface only reduces greatly the memory footprint (sort of like the node-vertex structures in geometry). You mention there are options to keep this versatility, how difficult do you think would be to implement them?

[pcarruscag]

Hi @rsanfer,
I was afraid you were going to say that...
On the subject of memory savings, and taking the specific case of CFEAVariable, the class is 304 bytes before it allocates data, once data is allocated for a plain nonlinear 3D problem the effective size is 704 bytes, looking at the constructors, for such a problem we allocate 22 doubles or 176 bytes. So effectively 75% of the memory (equivalent to 66 doubles) is wasted.
It is important to note that this "efficiency" figure does not get much better even when this class allocates all it members (for which one needs fsi + adjoint + unsteady). This happens because operator new [] needs to store the size of the allocation together with the data for when delete [] is called so you get "pointer + size + data" rounded to a multiple of 16B because of the default alignment.
Therefore my opinion is that we can "over allocate" and still save memory, unless you need massive amounts of data for the boundary points(?)... But looking at CFEABoundVariable I don't think that is the case.
If we want to be optimum in terms of performance and memory usage, we would need to mimic the CGeometry strategy in CSolver.
If we do not want to mimic that structure we need to either be able to rely on some ordering of the special cases, e.g. boundary variables stored first, so that at runtime:

double* Get(int iNode) {
    if (iNode < nBoundary)
        return &special[iNode];
    else
        return &normal[iNode-nBoundary];
}

If we cannot rely on some order we would need to tag nodes "special" or "normal" and then use some form of map for special cases... All this would cost a lot of indirection.

[rsanfer]

Hi @pcarruscag, thanks for the quick response.

There are still several possible memory reductions on the FEA solver by moving stuff to the vertices. Even though over-allocating would probably be fine for the structural solver as problems remain relatively small, I am replicating this structure for the mesh solver on the fluid side, and in there we should be careful not to add tons of doubles that are only necessary at the boundaries to the whole variable structure. I might also need a similar thing in the flow solver, e.g., to store boundary loads for the adjoint solver.

I believe that relying on some ordering would be quite complex and require a lot of reworking of the code, everything in the solvers links directly to the point index, which gets renumbered at the beginning. So forcing the vertex nodes to go on to the initial positions doesn't seem like an easy one. And, as you say, maybe if we flag the nodes (as boundary or interior) we end up pretty much at the starting point.

Possibly, the best alternative would be to replicate the structure in CGeometry in CSolver, therefore adding some CVertexVariable parent class and child for the different solvers. There are other boundary magnitudes that could be interesting to store as well. It would require some work, but given that there is a direct correspondence iMarker+iVertex -> iPoint, I would expect that to be the least invasive approach.

[pcarruscag]

Ok, so maybe we can start with over-allocating, and then move to a CVertexVariable approach?
If the extra allocations are for a few tractions/displacements we will save memory anyway.

[pcarruscag]

On the subject of virtual functions I would like to put another idea forward.
After moving things around in #725 I noticed that we have tons of inline virtual methods.
The keyword inline has two meanings to the compiler:

• "Dear merciful compiler please copy paste the body of this function and then do all your wonderful optimizations, if that pleases your excellency."
• "Dear forgiving compiler, you will find this method defined in multiple units, please don't be mad" (i.e. ignore the one-definition-rule).

virtual means determine what version of the method to call at runtime. This is not compatible with the first (and often the intended one) meaning of inline, therefore the compiler will in general not inline those methods.
They will only be inlined if they are being called on a pointer to the derived class that does not declare the method to be virtual anymore. CSolver knows what variables it creates and so in hot areas of the code it could do a static downcast to allow inlining (e.g. static_cast<CEulerVariable*>(node[iNode])->DoStuff()).
Where is this important? For example when computing gradients, where simple additions and subtractions are hidden behind virtual functions.
If you are worried about maintenance each solver can typedef its most safe downcast level or better yet (or just more modern), methods that could benefit from this can be templated for the type of downcast.
Those in favour say Yea those against say Nay.

[economon]

Aok on inlines.

Looking forward to the memory layout concept.. have you also thought about how we might combine this with vectorization (SIMD) across the edges in the residual / gradient / limiter routines? There's probably another 2x speedup available for those kernels.

Lastly, I don't think we should throw out the possibility of removing the CVariable class entirely in favor of allocating memory directly in the solver classes. We have thought about this in the past, and it is not used deliberately in the DG solver, for example. Food for thought.

[pcarruscag]

I've been doing some recreational reading on strategies that people have used for vectorization of fluxes but so far I have not thought of a way to reconcile them with the nice encapsulation we have of the numerics. For gradients (at least GG) and limiters on the other hand I do have ideas. By the way is there a branch on the repo for the work you, Francisco, Intel, and others did on this kind of stuff?

Yeah making solution variables members of the solver would avoid a lot of virtual calls but it would make it harder to eventually reach the solution for post processing, or in other areas like Integration/Iteration where the solution is manipulated directly by non-solver classes, we will see...

[vdweide]

My experience is that nice looking code and performance are two conflicting items, unfortunately. Apart from being contiguous in memory, you only obtain good performance when the vector length is a multiple of 8 or 16. This typically leads to padding the array, which does not improve readability. Furthermore, the best optimization is obtained when the vector length is known at compile time, typically leading to even uglier looking constructions,

[economon]

I can share that version of the code.. but beware, just as Edwin mentions, it is nasty looking and hard to compile (there are intrinsics and all sorts of architecture-specific changes). But, interesting to look at, of course. There you will see that we basically moved up all of the CVariable data (reordered contiguously) and methods to the solver to eliminate the extra layer of indirection, while still keeping access for outside classes, i.e. something like solver[FLOW_SOL]->node[iPoint]->GetPressure() becomes solver[FLOW_SOL]->GetPressure(iPoint).

Since you mentioned it, another target is the CNumerics classes. The flexibility is nice, but does it really need an entire set of classes with getters/setters/virtual functions? I think that most of the ComputeResidual() functions could be moved up as methods in the solver classes as well to make more efficient residual kernel loops.

In the early days, we thought completely removing CVariable and CNumerics would be very intrusive to change and would hurt flexibility. While the former is still true (although certainly doable), the latter does not seem to me as high of a priority any longer. Over the years, I don't think that folks have really modified the CNumerics classes all that often (once you have implemented the Roe method, it is more or less done). In addition, the CVariable and CNumerics child classes all live almost entirely within their own solvers, so we aren't taking much advantage of code reuse to justify having them as separate classes. For instance, if we just copied the relevant ComputeResidual() routines up into their solvers from numerics, we would have very little code duplication in the end, which is an indicator that we may not really need the extra baggage.

[clarkpede]

In the early days, we thought completely removing CVariable and CNumerics would be very intrusive to change and would hurt flexibility... Over the years, I don't think that folks have really modified the CNumerics classes all that often (once you have implemented the Roe method, it is more or less done).

As someone who has worked with the CNumerics classes a lot, I can say that we really don't gain much by keeping them as separate classes. On the surface, it seems like a good use of the strategy pattern. But from a pragmatic approach, it's not as flexible as OO purists would like you to believe. First off, if you need a new variable in the CNumerics class, you need to change both the solver class and the numerics class, since the solver class calls Setters to set the numerics variables. If you don't call the right Setters, you get a silent failure (unless you're using valgrind). The passing of variables between solver and numerics also leads to a lot of code duplication. And for many solvers (e.g. the viscous NS numerics or the source terms in the SST model) you only end up with one numerics class anyways.

I also imagine there's a strong performance hit from all the Setters and copying of variables. But I haven't taken the time to check for myself.

In my opinion, you could keep 80% of the current flexibility by using different functions instead of different classes. But I also think we would have to work carefully, or else we'll make the code harder to maintain.