Large-data memory layout for equation systems

Data layout requirements and design

How should field values, storing particle properties (in Walker) should be stored in memory? How should mesh field data (associated to nodes, elements, faces, etc., in Inciter) should be stored in memory? These are the single largest chunks of data a particle-, and/or mesh-based code operates on. The data layout, i.e., how the data is stored and organized in memory, determines how the data is accessed and potentially has a first-degree effect on overall performance.

Possibilities

1. Unknown-major, in which various physical properties, e.g., position, velocity, energy, etc., i.e., the unknowns in a solver, of a single particle or mesh entity are close to each other in memory. For example:

   where the are governed by one equation (e.g., position), the are governed by another equation (e.g., velocity), and the are governed by a third equation (e.g., energy), etc. Here the first letter denotes a physical quantity, while the second is the particle number or mesh entity index. If the algorithm advances the properties at the unknowns one equation at a time, data in memory is accessed by having to jump a distance that corresponds to the number of scalar physical variables per particle or mesh entity. In the example, the update will have to jump as are updated.

2. Property-major, in which the same type of physical properties are close to each other in memory. For example,

   The legend here is the same as in unknown-major: the first letter denotes a physical quantity, while the second is the particle or mesh entity index. If the algorithm advances the properties at the unknowns one equation at a time, data in memory is accessed contiguously in memory as the properties are contiguously stored.

Discussion

A property-major storage, case 2 above, seems to be the most efficient at first sight, as it stores data, as it is read and written by the equation algorithms, contiguously. However, data access is contiguous only if the particle properties (or data stored at mesh entities) are independent, i.e., if there is no coupling among the equations. Unfortunately, this is rarely the case, at least not for fluid dynamics or any computational physics problem that is interesting enough. For example in a Lagrangian particle code, velocity is used by the position update, and velocity is required by the energy update. The same is true for a mesh-based solver, where the physical variables are coupled, i.e., their update needs other physical variables at the same mesh entity (and usually that of the neighbors). Depending on the physical approximation, density (or mass) may be required for all equations. The stronger the equations are coupled, the more far-reads are required for a given update with a property-major data layout. These far-reads are practically always cache misses, as the property-major storage stores the physical variables for the same particle or mesh entity very far in memory, e.g., the distance between and is the number of particles or mesh entities. While the unknown-major storage, case 1 above, inherently stores data non-contiguously, the distance between properties of a single particle (or a single mesh entity) is relatively small, i.e., the number of properties, which may incur less cache misses as several particles (or nodes) but all of their physical properties at the same unknown could fit into cache.

Assuming strong coupling among the variables, the unknown-major storage will be favored, but it would be nice if the design allowed for both layouts, so depending on the type of equations the most appropriate layout could be selected. If such a design is maintainable, there is still a question whether the data layout selection should be done at compile-, or run-time.

Assessment of the Blaze library that offers a similar choice

Have looked at https://bitbucket.org/blaze-lib/blaze which implements row-, and column-major matrix classes based on a template argument. See, e.g., <blaze-1.5>/blaze/math/dense/StaticMatrix.h, which reveals that the template argument (bool) SO selects between row-, or column-major internal storage. Then SO is used at both compile-time (e.g., by the class-user, when instantiating the type of the matrix), as well as run-time (e.g., the implementation of isDefault()). Both compile-time and run-time usage of the SO template arguments are problematic:

• The compile-time usage duplicates a lot of code by having to provide similar implementations for the element-access operator() of StaticMatrix specialized to column-major. There is a generic implementation for SO for everything that is agnostic of SO, and there is a specialization when SO is column-major.
• The run-time usage also duplicates code by doing an if-test on SO in, e.g., isDefault(). Is there a better way of doing this? If there are only two types of data layout (unknown-, and property-major), code duplication should not be too much of an issue. However, the implementation of unknown-property data access must be absolutely zero run-time cost. This means the selection must be at compile-time and the element access must be absolutely invisible to client code. In other words, there must be no re-implementation of a time-integrator for an equation just because the data access is different.

Assessment of the Kokkos library that offers a similar choice

Since the first implementation of the configurable data layout, described below, the Kokkos library, https://github.com/kokkos/kokkos , from Sandia National Labs, has been released. Kokkos appears to have all the requirements described here, and many more other features. It does have compile-time configurable data layouts, views, for the purpose of optimal data access on various compute devices, such as multi-core CPUs, many-core accelerators, and Graphics Processing Units. Kokkos appears to provide an excellent abstraction for data layout abstraction. However, Kokkos provides a lot more functionality than just data layout abstraction, e.g., it can generate low level code for various devices using its configurable views. The currently available Kokkos back-ends are OpenMP, pthreads, and CUDA. Thus Kokkos provides abstractions for shared-memory parallelism. While Kokkos has been successfully used in Charm++ code, at this point (Jan 2016) we opt for not adopting Kokkos' views and keep our original data layout abstractions, discussed below. The reasons:

• At this point we are too early in the development of Quinoa and its tools to be able to say that the current abstractions Charm++ provides are sufficient or not to strike the correct balance between performance and productivity. In particular, it could be possible that Charm++'s abstractions (using overdecomposition), which we use already, are enough. In that case, using a single abstraction for parallelism is preferable to two.
• At this point we would really only need the compile-time configurable data layout abstractions from Kokkos and not the other rich features, such as abstractions for loop-level shared-memory parallelism.

As far I can tell, the views of an array in Kokkos provide a very similar abstraction for data layout and thus memory access what is described below.

Requirements

Zero-cost is achieved via type-based compile-time polymorphism. This is controlled via a cmake variable.

The particle (or mesh entity) data is a logically 3-dimensional array that stores the particle properties or unknowns at mesh entities, e.g., faces, cell centers, nodes, etc.

For clarity, the discussion below will use the expression "particle" for a particle-code and will not specifically mention mesh-based unknowns, stored at the nodes, elements, faces, etc., of a computational mesh. The former is for a particle-code, the latter is for a mesh-code. The data layout discussion is independent of whether particles or mesh entities (nodes, cells, etc.) are used, and accordingly both Walker and Inciter uses the same implementation.

In principle there are a total of 6 permutations:

1. ParEqComp: [ particle ] [ equation ] [ component ]
2. ParCompEq: [ particle ] [ component ] [ equation ]
3. EqCompPar: [ equation ] [ component ] [ particle ]
4. EqParComp: [ equation ] [ particle ] [ component ]
5. CompEqPar: [ component ] [ equation ] [ particle ]
6. CompParEq: [ component ] [ particle ] [ equation ]

Of these 6 we only consider those where component follows equation. (For those layouts where equation follows component the access would be unnecessarily complicated by the potentially unequal number of components for different equations which is not known at compile-time and thus does not allow some optimizations.) This decision leaves us with the following choices:

1. ParEqComp: [ particle ] [ equation ] [ component ]
3. EqCompPar: [ equation ] [ component ] [ particle ]
4. EqParComp: [ equation ] [ particle ] [ component ]

Access is based on the 3 coordinates: particle (or unknown), component, and offset. Particle is the particle index, component denotes the given component of a vector equation, e.g., velocity has 3 components, a multi-material turbulent mix model governed by the Dirichlet SDE has scalars (components), and offset is determined by the relative position of the given equation compared to the other equations. Using these 3 coordinates the index calculations for the above 3 cases are:

1. ParEqComp: [ particle ] [ equation ] [ component ]

     baseptr + particle*nprop + offset + component

where nprop is the total number of particle properties, e.g., 3 positions, 3 velocities, 5 scalars nprop = 11.

3. EqCompPar: [ equation ] [ component ] [ particle ]

     baseptr + (offset+component)*npar + particle

where npar is the total number of particles.

4. EqParComp: [ equation ] [ particle ] [ component ]

     baseptr + offset*npar + nce*particle + component

where nce is the number of components for the given equation. Since this would require another function argument (besides particle, component, and offset), and it costs an integer-multiply more than the other two layouts, we dismiss this layout, and only implement the following two:

1. ParEqComp - Particle-major
3. EqCompPar - Equation-major

These options are exposed via a cmake variable and can be switched before a build.

Data layout implementation and assembly

This section documents the implementation and the assembly code, produced by the compilers, of the compile-time configurable data-access policy discussed above. The implementation is via a thin data-access interface with zero run-time cost, no code-duplication, and in a way that is invisible to client code.

const bool ParticleMajor = true;
const bool PropertyMajor = false;

template< bool Major >
class Data {

  private:
   // Transform a compile-time bool into a type
   template< bool m >
   struct int2type {
     enum { value = m };
   };

   // Overloads for particle-, and property-major accesses
   inline
   tk::real& access( int particle, int property, int2type<ParticleMajor> ) {
     return *(m_ptr + particle*m_nprop + m_offset + property);
   }
   inline
   tk::real& access( int particle, int property, int2type<PropertyMajor> ) {
     // This is the same for now, not called, irrelevant in zero-cost-test
     return *(m_ptr + particle*m_nprop + m_offset + property);
   }

   tk::real* const m_ptr;
   const int m_nprop;
   const int m_offset;

  public:
    // Constructor
    Data( tk::real* const ptr, int nprop, int offset ) :
      m_ptr(ptr), m_nprop(nprop), m_offset(offset) {}

    // Access dispatch
    inline tk::real& operator()( int particle, int property ) {
      return access( particle, property, int2type<Major>() );
    }
};

Data< Layout > d( particles, m_nprop, m_offset );
Model::aa = d( 34, 3 );
Model::bb = *(m_particles + 34*m_nprop + m_offset + 3);

Model {
  ...
  public:
    tk::real aa;
    tk::real bb;
  ...
}

std::cout << m_mix->aa << m_mix->bb;

     .loc    143 42 20  ; <quinoa>/src/DiffEq/Dirichlet.h:42:20
.Ltmp27038:
     lea     RDI, QWORD PTR [RBP - 56]
     mov     ESI, 34
     mov     EDX, 3
     call    _ZN6quinoa4DataILb1EEclEii
.Ltmp27039:
     mov     QWORD PTR [RBP - 176], RAX # 8-byte Spill
     jmp     .LBB2550_7
.LBB2550_7:
     mov     RAX, QWORD PTR [RBP - 176] # 8-byte Reload
     movsd   XMM0, QWORD PTR [RAX]
     mov     RCX, QWORD PTR [RBP - 64] # 8-byte Reload
     movsd   QWORD PTR [RCX + 8], XMM0
     .loc    143 43 0    ; <quinoa>/src/DiffEq/Dirichlet.h:43:0
     mov     RDX, QWORD PTR [RCX + 32]
     imul    ESI, DWORD PTR [RCX + 48], 34
     movsxd  RDI, ESI
     shl     RDI, 3
     add     RDX, RDI
     movsxd  RDI, DWORD PTR [RCX + 52]
     movsd   XMM0, QWORD PTR [RDX + 8*RDI + 24]
     movsd   QWORD PTR [RCX + 16], XMM0

.loc    144 42 20  ; <quinoa>/src/DiffEq/Dirichlet.h:42:20
movsd   XMM0, QWORD PTR [R14 + 8*RAX + 24]
movsd   QWORD PTR [R13 + 8], XMM0
.loc    144 43 0    ; <quinoa>/src/DiffEq/Dirichlet.h:43:0
mov     RCX, QWORD PTR [R13 + 32]
movsd   XMM0, QWORD PTR [RCX + 8*RAX + 24]
movsd   QWORD PTR [R13 + 16], XMM0

Data layout benchmark

This section documents the benchmark of the implementation of the compile-time configurable data-access policy discussed above. The implementation is via a thin data-access interface with zero run-time cost, no code-duplication, and in a way that is invisible to client code.

walker

  nstep 100000  # Max number of time steps
  term  140.0   # Max time
  dt    0.05    # Time step size
  npar  40000   # Number of particles

  ttyi  100     # TTY output interval

  rngs
    mkl_mrg32k3a seed 0 end
  end

  dirichlet
    ncomp 100  # = K = N-1
    b     0.1 1.5 0.1 1.5 0.1 1.5 0.1 1.5 0.1 1.5
          0.1 1.5 0.1 1.5 0.1 1.5 0.1 1.5 0.1 1.5
          0.1 1.5 0.1 1.5 0.1 1.5 0.1 1.5 0.1 1.5
          0.1 1.5 0.1 1.5 0.1 1.5 0.1 1.5 0.1 1.5
          0.1 1.5 0.1 1.5 0.1 1.5 0.1 1.5 0.1 1.5
          0.1 1.5 0.1 1.5 0.1 1.5 0.1 1.5 0.1 1.5
          0.1 1.5 0.1 1.5 0.1 1.5 0.1 1.5 0.1 1.5
          0.1 1.5 0.1 1.5 0.1 1.5 0.1 1.5 0.1 1.5
          0.1 1.5 0.1 1.5 0.1 1.5 0.1 1.5 0.1 1.5
          0.1 1.5 0.1 1.5 0.1 1.5 0.1 1.5 0.1 1.5
    end
    S     0.625 0.4 0.625 0.4 0.625 0.4 0.625 0.4 0.625 0.4
          0.625 0.4 0.625 0.4 0.625 0.4 0.625 0.4 0.625 0.4
          0.625 0.4 0.625 0.4 0.625 0.4 0.625 0.4 0.625 0.4
          0.625 0.4 0.625 0.4 0.625 0.4 0.625 0.4 0.625 0.4
          0.625 0.4 0.625 0.4 0.625 0.4 0.625 0.4 0.625 0.4
          0.625 0.4 0.625 0.4 0.625 0.4 0.625 0.4 0.625 0.4
          0.625 0.4 0.625 0.4 0.625 0.4 0.625 0.4 0.625 0.4
          0.625 0.4 0.625 0.4 0.625 0.4 0.625 0.4 0.625 0.4
          0.625 0.4 0.625 0.4 0.625 0.4 0.625 0.4 0.625 0.4
          0.625 0.4 0.625 0.4 0.625 0.4 0.625 0.4 0.625 0.4
    end
    kappa 0.0125 0.3 0.0125 0.3 0.0125 0.3 0.0125 0.3 0.0125 0.3
          0.0125 0.3 0.0125 0.3 0.0125 0.3 0.0125 0.3 0.0125 0.3
          0.0125 0.3 0.0125 0.3 0.0125 0.3 0.0125 0.3 0.0125 0.3
          0.0125 0.3 0.0125 0.3 0.0125 0.3 0.0125 0.3 0.0125 0.3
          0.0125 0.3 0.0125 0.3 0.0125 0.3 0.0125 0.3 0.0125 0.3
          0.0125 0.3 0.0125 0.3 0.0125 0.3 0.0125 0.3 0.0125 0.3
          0.0125 0.3 0.0125 0.3 0.0125 0.3 0.0125 0.3 0.0125 0.3
          0.0125 0.3 0.0125 0.3 0.0125 0.3 0.0125 0.3 0.0125 0.3
          0.0125 0.3 0.0125 0.3 0.0125 0.3 0.0125 0.3 0.0125 0.3
          0.0125 0.3 0.0125 0.3 0.0125 0.3 0.0125 0.3 0.0125 0.3
    end
    rng mkl_mrg32k3a
  end

  statistics    # Estimate statistics
    <Y1>        # mean of Y1
    <Y2>
    <y1y1>      # variance of Y1 = <(Y1-<Y1>)^2> = <y1^2>
    <y2y2>
    <y1y2>
  end

end

//! Advance particles
void advance(int p, int tid, tk::real dt) override {
  // Get access to particle scalars
  tk::real* y = m_particles.ptr() + p*m_nprop + m_offset;
  
  // Compute Nth scalar
  tk::real yn = 1.0 - y[0];
  for (int i=1; i<m_ncomp; ++i) yn -= y[i];
  
  // Generate Gaussian random numbers with zero mean and unit variance
  tk::real dW[m_ncomp];
  m_rng->gaussian( tid, m_ncomp, dW );
  
  // Advance first m_ncomp (K=N-1) scalars
  for (int i=0; i<m_ncomp; ++i) {
    tk::real d = m_k[i]*y[i]*yn*dt;
    if (d > 0.0) d = sqrt(d); else d = 0.0;
    y[i] += 0.5*m_b[i]*(m_S[i]*yn - (1.0-m_S[i])*y[i])*dt + d*dW[i];
  }
}

//! Advance particles
void advance(int p, int tid, tk::real dt) override {
  // Compute Nth scalar
  tk::real yn = 1.0 - m_particles(p, 0, m_offset);
  for (int i=1; i<m_ncomp; ++i) yn -= m_particles(p, i, m_offset);

  // Generate Gaussian random numbers with zero mean and unit variance
  tk::real dW[m_ncomp];
  m_rng->gaussian( tid, m_ncomp, dW );

  // Advance first m_ncomp (K=N-1) scalars
  for (int i=0; i<m_ncomp; ++i) {
    tk::real d = m_k[i] * m_particles(p, i, m_offset) * yn * dt;
    if (d > 0.0) d = sqrt(d); else d = 0.0;
    m_particles(p, i, m_offset) +=
      0.5*m_b[i]*(m_S[i]*yn - (1.0-m_S[i]) * m_particles(p, i, m_offset) )*dt
      + d*dW[i];
  }
}

Conclusion

As the implementation is not a significant performance hit, the equation advancement algorithms and particle-, and mesh-property data access are used only via the data-layout-independent interface. The data layout can be changed at compile time. Data access is abstracted (and optimized) away.

Note that the above discussion is independent of whether a particle-code or a mesh-code is used. From the data layout viewpoint the particle-major and mesh-field-major are equivalent and in the code, and in tk::Data, this is called the unknown-major. See src/Base/Data.h and its specializations in src/Base/Fields.h and src/Base/Particles.h. The data layouts for mesh-fields and particles can be configured independently at compile time via cmake variables PARTICLE_DATA_LAYOUT and FIELD_DATA_LAYOUT, see also cmake/ConfigureDataLayout.cmake.

For the API, see tk::Data, and for the full (and current) implementation, see src/Base/Data.h.