Creating a SOA with nice OOP like interface

[Zachary 42!]

Hi everyone,

I saw on lecture 22 “Some data structure design considerations” Prof. Wolfgang described a struct of arrays (SOA) with an accessor interface.  This keeps nice encapsulation while also providing nice cache use.  I would like to understand this design better for my personal programming growth, separating interface and data so the data is stored optimally.  

Here is a snippet from the iterators section of the documentation:

"Since dereferencing iterators yields accessor objects, these calls are to member functions Accessor::vertex(), Accessor::child() etc. These in turn figure out the relevant data from the various data structures that store this data. How this is actually done and what data structures are used is not really of concern to authors of applications in deal.II. In particular, by hiding the actual data structures we are able to store data in an efficient way, not necessarily in a way that makes it easily accessible or understandable to application writers.”

What classes should I look at that shows this design the most succinctly?  How is the data actually stored?  This design seems very important for writing data oriented code but still provided a nice interface layer on top and I would like to get comfortable with it.

Cheers,

Zachary

[peterrum]

Hi Zachary,

I haven't watched the lecture, but I think I can, nevertheless, answer this question since I have been refactoring and generalizing the relevant data structures in the last months. 

The key aspect is that the data is separated from the class that gives coordinated access to the data. In our case, the `Triangulation` classes (and similarity the `DoFHandler`) return iterators, which are more or less wrapper around `TriaAccessor` or `CellAccessor` (which is also a `TriaAccessor` with some additional functionalities). The core of these two classes are the internal fields `present_level` and `present_index`, which function as pointers into some arrays (within the `Triangulation` class).

Let's for example have a look at the implementation of `CellAccessor::subdomain_id()`:

```cpp

template <int dim, int spacedim>
inline types::subdomain_id
CellAccessor<dim, spacedim>::subdomain_id() const
{
  return this->tria->levels[this->present_level]->subdomain_ids[this->present_index];
}

```

which access via the `Triangulation` class `internal::TriangulationImplementation::TriaLevel::subdomain_ids` array. The functions `child()` and `vertex()` are setup with a similar logic but are a bit more complicated: the first one has a stride of 2^dim and the second one uses many indirections (e.g. in the case of 3D cell, querying the vertex of a face, which queries the vertex of a line - however each query is similar to the one of `subdomain_id()`).

The assumption of this storage is that in one place of your code you loop over all cells and only need one or a few properties of a cell (e.g. `subdomain_id` or `boundary_id`). As a consequence, how to store data is very much determined how you data will be processed. 

If you actually want to look at the internal data structures of the `Triangulation` class with the relevant SOA structures `TriaLevel`, `TriaFaces`, TriaObjects` (and the accessor classes), the following UML diagram might be of help: https://github.com/dealii/dealii/issues/10582#issue-644048824

Peter

[Zachary Streeter]

Hi Peter,

Thank you for your email.  Here is what I have so far:

```c++

// This will be filled in by looping over number of orbitals, with each orbital containing an array of coefs.

struct Orbital_Data{

\\TODO: need a constructor that allocates memory for all the orbital's info here?

private:

  std::shared_ptr<std::complex<double>[ ] > orbital_coefs;

  std::shared_ptr<uint42_t[ ]> orbital_coefs_index;

}

struct Orbital{

\\TODO: How to handle the ith_orbital index?  

\\ Something like: Orbital(const uint32_t &a_ith_orbital)

  inline std:complex<double>

  coef(uint32_t i) const {

    uint32_t coef_index =

      orbital_data->orbital_coefs_index[ith_orbital*size_of_coefs_array + i]; // Note the ith_orbital here used as a stride

  return orbital_data->orbital_coefs[ coefs_index];

  }

private:

  Orbital_Data *orbital_data;

  uint32_t ith_orbital;  // This is the heart of SOA

}

```

I would rather allocate enough space instead of using std::vectors in a constructor possible.  Does this look reasonable to you?

My two questions are my TODOs.  I'm thinking about passing in the size the arrays need to be into a constructor for Orbital_Data.  I am looking to fill in this data inside a loop over orbitals.  So I need to know the ith_orbital and pass it into the Orbital struct somehow.

This follows Prof. Wolfgang's example in his lecture.  I am just missing a few little things conceptually.

Cheers,

Zachary

[peterrum]

Hi Zachary,

Regarding your two TODOs. I am explaining how we do it in deal.II and then you can decide if this approach also fits your needs.

Index handling

The simple case is that the indices (level, index within the level) are passed to `TriaAccessorBase` in the constructor (like here https://github.com/dealii/dealii/blob/a56c8585863b774de9dce5f7dd5c334ef56f3a51/include/deal.II/grid/tria_accessor.h#L341-L344). The more common approach is that `TriaAccessorBase` is used in the context of an iterator; the iterator calls the `TriaAccessorBase::operator++()` function, which increments the indices internally appropriately (see: https://github.com/dealii/dealii/blob/a56c8585863b774de9dce5f7dd5c334ef56f3a51/include/deal.II/grid/tria_accessor.templates.h#L198-L236).

Allocation

In deal.II we have dedicated allocation functions, called `reserve_space` see:  https://github.com/dealii/dealii/blob/a56c8585863b774de9dce5f7dd5c334ef56f3a51/source/grid/tria.cc#L2638-L2744. But doing the reservation in the constructor is also a way if you now the size at time of construction, which we don't know in the Triangulation class, and it does not change over time (e.g., due to adaptivity).

A general remark: it might be also useful - depending on the regularity of data-access pattern -  to have a look into explicit SIMD vectorization, i.e., in your case to pack 8 orbitals in the case of AVX512 together. For this purpose, deal.II has the VectorizedArray struct, which is a wrapper class around intrinsic instruction so that one does not have mess with these.

I hope this helps. 

Peter