CD_ParticleSoA.H

Go to the documentation of this file.

/*
 * SPDX-FileCopyrightText: 2021-2026 SINTEF Energy Research
 *
 * SPDX-License-Identifier: GPL-3.0-or-later
 */
 
#ifndef CD_PARTICLESOA_H
#define CD_PARTICLESOA_H
 
// Std includes
#include <algorithm>
#include <array>
#include <cstddef>
#include <cstdint>
#include <cstdlib>
#include <cstring>
#include <limits>
#include <tuple>
#include <type_traits>
#include <utility>
#include <vector>
 
// Chombo includes (kept minimal -- only the grid-interop geometry types; storage uses
// fundamental C++ types so the container is not tied to Chombo's Real precision).
#include <SPACE.H>
#include <RealVect.H>
#include <IntVect.H>
#include <Box.H>
 
// Our includes
#include <CD_ParticleMemory.H>
#include <CD_NamespaceHeader.H>
 
#ifndef CD_PARTICLE_REAL
#define CD_PARTICLE_REAL double
#endif
 
using ParticleReal = CD_PARTICLE_REAL;
 
using ParticleID = std::int64_t;
 
using RankID = std::int32_t;
 
enum class GhostType : std::uint8_t
{
  Valid     = 0, 
  SameLevel = 1, 
  Coarse    = 2, 
  Fine      = 3  
};
 
template <typename P>
struct ParticleTraits;
 
namespace H5Part {
 
template <typename M>
struct Scalar
{
  static constexpr bool s_isVector = false;
 
  const char* name;
 
  M member;
};
 
template <typename M>
Scalar(const char*, M) -> Scalar<M>;
 
template <typename... Ms>
struct Vector
{
  static constexpr bool s_isVector = true;
 
  const char* name;
 
  std::tuple<Ms...> members;
 
  constexpr Vector(const char* a_name, Ms... a_members) noexcept : name(a_name), members(a_members...)
  {}
};
 
template <typename... Ms>
Vector(const char*, Ms...) -> Vector<Ms...>;
 
} // namespace H5Part
 
struct NoPayload
{
};
 
template <>
struct ParticleTraits<NoPayload>
{
  static constexpr auto columns = std::make_tuple();
};
 
namespace detail {
 
template <typename M>
struct MemberPointerTraits;
 
template <typename C, typename T>
struct MemberPointerTraits<T C::*>
{
  using ClassType = C;
  using ValueType = T;
};
 
template <typename A, typename B>
constexpr bool
eqMember(A, B) noexcept
{
  return false;
}
 
template <typename A>
constexpr bool
eqMember(A a_lhs, A a_rhs) noexcept
{
  return a_lhs == a_rhs;
}
 
template <typename Tup, typename M>
constexpr std::size_t
indexOfImpl(const Tup&, M, std::index_sequence<>) noexcept
{
  return std::tuple_size<Tup>::value;
}
 
template <typename Tup, typename M, std::size_t I, std::size_t... Rest>
constexpr std::size_t
indexOfImpl(const Tup& a_t, M a_target, std::index_sequence<I, Rest...>) noexcept
{
  return eqMember(std::get<I>(a_t), a_target) ? I : indexOfImpl(a_t, a_target, std::index_sequence<Rest...>{});
}
 
template <typename Tup, typename M>
constexpr std::size_t
indexOf(const Tup& a_t, M a_target) noexcept
{
  return indexOfImpl(a_t, a_target, std::make_index_sequence<std::tuple_size<Tup>::value>{});
}
 
template <typename T, std::size_t N>
struct RepeatTuple
{
  template <std::size_t>
  using Always = T;
 
  template <std::size_t... I>
  static std::tuple<Always<I>...> make(std::index_sequence<I...>);
 
  using Type = decltype(make(std::make_index_sequence<N>{}));
};
 
template <typename Tuple>
struct ValueTuple;
 
template <typename... Ms>
struct ValueTuple<std::tuple<Ms...>>
{
  using Type = std::tuple<typename MemberPointerTraits<Ms>::ValueType...>;
};
 
constexpr std::size_t
sumSizes() noexcept
{
  return 0;
}
 
template <typename... Rest>
constexpr std::size_t
sumSizes(std::size_t a_first, Rest... a_rest) noexcept
{
  return a_first + sumSizes(a_rest...);
}
 
template <typename Tuple>
struct SumSizeof;
 
template <typename... Ts>
struct SumSizeof<std::tuple<Ts...>>
{
  static constexpr std::size_t s_value = sumSizes(sizeof(Ts)...);
};
 
template <typename Tuple>
struct AllTriviallyCopyable;
 
template <typename... Ts>
struct AllTriviallyCopyable<std::tuple<Ts...>>
{
  static constexpr bool s_value = (std::is_trivially_copyable<Ts>::value && ...);
};
 
template <typename...>
struct MakeVoid
{
  using Type = void;
};
 
template <typename... Ts>
using Void = typename MakeVoid<Ts...>::Type;
 
template <typename Traits, std::size_t... I>
constexpr std::index_sequence<indexOf(Traits::columns, std::get<I>(Traits::h5Columns))...>
deriveH5Seq([[maybe_unused]] std::index_sequence<I...> a_seq) noexcept
{
  return {};
}
 
template <typename Traits, typename Fallback, typename = void>
struct H5ColumnSubset
{
  using Type = Fallback;
};
 
template <typename Traits, typename Fallback>
struct H5ColumnSubset<Traits, Fallback, Void<decltype(Traits::h5Columns)>>
{
  using H5Tuple = std::remove_cv_t<std::remove_reference_t<decltype(Traits::h5Columns)>>;
  using Type = decltype(deriveH5Seq<Traits>(std::make_index_sequence<std::tuple_size<H5Tuple>::value>{}));
};
 
template <typename T>
inline void
pushBytes(unsigned char*& a_buf, const T& a_value) noexcept
{
  static_assert(std::is_trivially_copyable<T>::value, "SoA columns must be trivially copyable");
  std::memcpy(a_buf, &a_value, sizeof(T));
  a_buf += sizeof(T);
}
 
template <typename T>
inline void
pullBytes(const unsigned char*& a_buf, T& a_value) noexcept
{
  static_assert(std::is_trivially_copyable<T>::value, "SoA columns must be trivially copyable");
  std::memcpy(&a_value, a_buf, sizeof(T));
  a_buf += sizeof(T);
}
 
} // namespace detail
 
template <typename P = NoPayload, typename Traits = ParticleTraits<P>>
class ParticleSoA
{
public:
 
  using PayloadColumns = std::remove_cv_t<std::remove_reference_t<decltype(Traits::columns)>>;
 
  static constexpr std::size_t s_numPayloadColumns = std::tuple_size<PayloadColumns>::value;
 
  static constexpr std::size_t s_weightCol = SpaceDim;
 
  static constexpr std::size_t s_idCol = SpaceDim + 1;
 
  static constexpr std::size_t s_rankCol = SpaceDim + 2;
 
  static constexpr std::size_t s_ghostCol = SpaceDim + 3;
 
  static constexpr std::size_t s_payloadBegin = SpaceDim + 4;
 
  static constexpr std::size_t s_numColumns = s_payloadBegin + s_numPayloadColumns;
 
  static constexpr std::size_t s_align = 64;
 
 
private:
  using PositionValues = typename detail::RepeatTuple<double, SpaceDim>::Type;
 
  using MandatoryValues = decltype(std::tuple_cat(std::declval<PositionValues>(),
                                                  std::declval<std::tuple<double, ParticleID, RankID, GhostType>>()));
 
  using PayloadValues = typename detail::ValueTuple<PayloadColumns>::Type;
 
public:
 
  using ColumnValueTuple = decltype(std::tuple_cat(std::declval<MandatoryValues>(), std::declval<PayloadValues>()));
 
  template <std::size_t K>
  using ColumnType = std::tuple_element_t<K, ColumnValueTuple>;
 
  template <typename M>
  static constexpr std::size_t
  columnIndex(const M a_member) noexcept
  {
    return s_payloadBegin + detail::indexOf(Traits::columns, a_member);
  }
 
  using H5PayloadColumns = typename detail::H5ColumnSubset<Traits, std::make_index_sequence<s_numPayloadColumns>>::Type;
 
 
  // ----- Sanity checks ------------------------------------------------------
  static_assert(std::is_default_constructible<P>::value, "payload type P must be default-constructible");
  static_assert(detail::AllTriviallyCopyable<ColumnValueTuple>::s_value,
                "every particle column must be trivially copyable for MPI/HDF5");
 
 
  ParticleSoA() = default;
 
  ~ParticleSoA()
  {
    ParticleMemory::removeBytes(ParticleMemory::Kind::Container, m_bytes);
 
    std::free(m_buffer);
  }
 
  ParticleSoA(const ParticleSoA&) = delete;
 
  ParticleSoA&
  operator=(const ParticleSoA&) = delete;
 
  ParticleSoA(ParticleSoA&& a_other) noexcept
  {
    this->steal(a_other);
  }
 
  ParticleSoA&
  operator=(ParticleSoA&& a_other) noexcept
  {
    if (this != &a_other) {
      ParticleMemory::removeBytes(ParticleMemory::Kind::Container, m_bytes);
 
      std::free(m_buffer);
      this->steal(a_other);
    }
    return *this;
  }
 
  [[nodiscard]] ParticleSoA
  deepCopy() const;
 
  void
  deepCopyTo(ParticleSoA& a_dst) const;
 
  void
  swap(ParticleSoA& a_other) noexcept
  {
    std::swap(m_buffer, a_other.m_buffer);
    std::swap(m_bytes, a_other.m_bytes);
    std::swap(m_capacity, a_other.m_capacity);
    std::swap(m_size, a_other.m_size);
    std::swap(m_sorted, a_other.m_sorted);
    std::swap(m_colPtrs, a_other.m_colPtrs);
    m_cellStart.swap(a_other.m_cellStart);
    std::swap(m_sortBox, a_other.m_sortBox);
    std::swap(m_sortDx, a_other.m_sortDx);
    std::swap(m_sortProbLo, a_other.m_sortProbLo);
  }
 
 
 
  std::size_t
  size() const noexcept
  {
    return m_size;
  }
 
  std::size_t
  capacity() const noexcept
  {
    return m_capacity;
  }
 
  bool
  empty() const noexcept
  {
    return m_size == 0;
  }
 
  void
  clear() noexcept
  {
    if (m_size > 0) {
      m_sorted = false;
    }
 
    m_size = 0;
  }
 
  void
  reserve(const std::size_t a_capacity);
 
  void
  resize(const std::size_t a_count);
 
  void
  shrinkToFit();
 
 
 
  void
  append(const RealVect& a_position, const double a_weight)
  {
    this->append(a_position, a_weight, P{});
  }
 
  void
  append(const RealVect& a_position, const double a_weight, const P& a_payload);
 
  void
  append(const ParticleSoA& a_other);
 
  void
  appendParticle(const ParticleSoA& a_src, const std::size_t a_index);
 
  void
  catenate(ParticleSoA& a_other)
  {
    if (this == &a_other || a_other.m_size == 0) {
      return;
    }
    if (m_size == 0) {
      this->swap(a_other); // zero-copy steal into an empty destination
      return;
    }
    this->append(a_other); // bulk per-column copy
    a_other.clear();       // move semantics: leave the source empty
  }
 
  P
  gather(const std::size_t a_index) const
  {
    CH_assert(a_index < m_size);
 
    return this->gatherPayload(a_index, std::make_index_sequence<s_numPayloadColumns>{});
  }
 
  void
  remove(const std::size_t a_index) noexcept
  {
    CH_assert(a_index < m_size);
 
    this->removeImpl(a_index, std::make_index_sequence<s_numColumns>{});
    m_size--;
    m_sorted = false;
  }
 
 
 
  template <std::size_t K>
  ColumnType<K>*
  columnByIndex() noexcept
  {
    return reinterpret_cast<ColumnType<K>*>(m_colPtrs[K]);
  }
 
  template <std::size_t K>
  const ColumnType<K>*
  columnByIndex() const noexcept
  {
    return reinterpret_cast<const ColumnType<K>*>(m_colPtrs[K]);
  }
 
  template <auto Member>
  auto*
  column() noexcept
  {
    return this->template columnByIndex<columnIndex(Member)>();
  }
 
  template <auto Member>
  const auto*
  column() const noexcept
  {
    return this->template columnByIndex<columnIndex(Member)>();
  }
 
  template <auto Member>
  auto&
  get(const std::size_t a_index) noexcept
  {
    CH_assert(a_index < m_size);
 
    return this->template columnByIndex<columnIndex(Member)>()[a_index];
  }
 
  template <auto Member>
  const auto&
  get(const std::size_t a_index) const noexcept
  {
    CH_assert(a_index < m_size);
 
    return this->template columnByIndex<columnIndex(Member)>()[a_index];
  }
 
  double*
  positionColumn(const int a_dir) noexcept
  {
    m_sorted = false;
 
    return reinterpret_cast<double*>(m_colPtrs[a_dir]);
  }
 
  const double*
  positionColumn(const int a_dir) const noexcept
  {
    return reinterpret_cast<const double*>(m_colPtrs[a_dir]);
  }
 
  double*
  weightColumn() noexcept
  {
    return this->template columnByIndex<s_weightCol>();
  }
 
  const double*
  weightColumn() const noexcept
  {
    return this->template columnByIndex<s_weightCol>();
  }
 
 
 
  RealVect
  position(const std::size_t a_index) const noexcept
  {
    CH_assert(a_index < m_size);
 
    RealVect x;
    for (int dir = 0; dir < SpaceDim; dir++) {
      x[dir] = this->positionColumn(dir)[a_index];
    }
    return x;
  }
 
  void
  setPosition(const std::size_t a_index, const RealVect& a_position) noexcept
  {
    CH_assert(a_index < m_size);
 
    // Writes through the non-const positionColumn(), which clears m_sorted.
    for (int dir = 0; dir < SpaceDim; dir++) {
      this->positionColumn(dir)[a_index] = a_position[dir];
    }
  }
 
  double&
  weight(const std::size_t a_index) noexcept
  {
    CH_assert(a_index < m_size);
 
    return this->template columnByIndex<s_weightCol>()[a_index];
  }
 
  const double&
  weight(const std::size_t a_index) const noexcept
  {
    CH_assert(a_index < m_size);
 
    return this->template columnByIndex<s_weightCol>()[a_index];
  }
 
  static constexpr ParticleID s_invalidID = -1;
 
  ParticleID&
  particleID(const std::size_t a_index) noexcept
  {
    CH_assert(a_index < m_size);
 
    return this->template columnByIndex<s_idCol>()[a_index];
  }
 
  const ParticleID&
  particleID(const std::size_t a_index) const noexcept
  {
    CH_assert(a_index < m_size);
 
    return this->template columnByIndex<s_idCol>()[a_index];
  }
 
  RankID&
  rankID(const std::size_t a_index) noexcept
  {
    CH_assert(a_index < m_size);
 
    return this->template columnByIndex<s_rankCol>()[a_index];
  }
 
  const RankID&
  rankID(const std::size_t a_index) const noexcept
  {
    CH_assert(a_index < m_size);
 
    return this->template columnByIndex<s_rankCol>()[a_index];
  }
 
  GhostType&
  ghost(const std::size_t a_index) noexcept
  {
    CH_assert(a_index < m_size);
 
    return this->template columnByIndex<s_ghostCol>()[a_index];
  }
 
  const GhostType&
  ghost(const std::size_t a_index) const noexcept
  {
    CH_assert(a_index < m_size);
 
    return this->template columnByIndex<s_ghostCol>()[a_index];
  }
 
  bool
  isGhost(const std::size_t a_index) const noexcept
  {
    return this->ghost(a_index) != GhostType::Valid;
  }
 
 
 
  static constexpr std::size_t
  bytesPerParticle() noexcept
  {
    return detail::SumSizeof<ColumnValueTuple>::s_value;
  }
 
  void
  linearizeParticle(void* a_buffer, const std::size_t a_index) const noexcept;
 
  void
  delinearizeAndAppend(const void* a_buffer);
 
 
 
  static constexpr std::size_t
  h5BytesPerParticle() noexcept
  {
    return (SpaceDim + 1) * sizeof(double) + h5PayloadSize(H5PayloadColumns{});
  }
 
  void
  h5LinearizeParticle(void* a_buffer, const std::size_t a_index) const noexcept;
 
  void
  h5DelinearizeAndAppend(const void* a_buffer);
 
 
 
  const void*
  data() const noexcept
  {
    return m_buffer;
  }
 
  void*
  data() noexcept
  {
    return m_buffer;
  }
 
  std::size_t
  byteSpan() const noexcept
  {
    if (m_size == 0) {
      return 0;
    }
    return this->byteSpanImpl(std::make_index_sequence<s_numColumns>{});
  }
 
 
 
  void
  sortByCell(const Box& a_box, const RealVect& a_dx, const RealVect& a_probLo);
 
  bool
  isSorted() const noexcept
  {
    return m_sorted;
  }
 
  bool
  isSortedAgainst(const Box& a_box, const RealVect& a_dx, const RealVect& a_probLo) const noexcept
  {
    return m_sorted && m_sortBox == a_box && m_sortDx == a_dx && m_sortProbLo == a_probLo;
  }
 
  std::size_t
  numCells() const noexcept
  {
    return m_cellStart.empty() ? 0 : m_cellStart.size() - 1;
  }
 
  std::size_t
  cellStart(const std::size_t a_cell) const noexcept
  {
    CH_assert(m_sorted);
    CH_assert(a_cell <= this->numCells());
 
    return m_cellStart[a_cell];
  }
 
  std::size_t
  particlesInCell(const std::size_t a_cell) const noexcept
  {
    CH_assert(m_sorted);
    CH_assert(a_cell < this->numCells());
 
    return m_cellStart[a_cell + 1] - m_cellStart[a_cell];
  }
 
  std::pair<std::size_t, std::size_t>
  cellRange(const std::size_t a_cell) const noexcept
  {
    CH_assert(m_sorted);
    CH_assert(a_cell < this->numCells());
 
    return {m_cellStart[a_cell], m_cellStart[a_cell + 1]};
  }
 
  void
  extractCell(const std::size_t a_cell, ParticleSoA& a_out) const
  {
    CH_assert(m_sorted);
 
    a_out.clear();
 
    const std::pair<std::size_t, std::size_t> range = this->cellRange(a_cell);
    for (std::size_t i = range.first; i < range.second; i++) {
      a_out.appendParticle(*this, i);
    }
  }
 
  const std::vector<std::size_t>&
  cellOffsets() const noexcept
  {
    CH_assert(m_sorted);
 
    return m_cellStart;
  }
 
  void
  adoptCellSort(const Box&               a_box,
                const RealVect&          a_dx,
                const RealVect&          a_probLo,
                std::vector<std::size_t> a_cellStart) noexcept
  {
    CH_assert(a_cellStart.size() == static_cast<std::size_t>(a_box.numPts()) + 1);
    CH_assert(a_cellStart.back() == m_size);
 
    m_cellStart  = std::move(a_cellStart);
    m_sortBox    = a_box;
    m_sortDx     = a_dx;
    m_sortProbLo = a_probLo;
    m_sorted     = true;
 
    CH_assert(this->checkCellSort(a_box, a_dx, a_probLo));
  }
 
 
protected:
  void* m_buffer = nullptr;
 
  std::size_t m_capacity = 0;
 
  std::size_t m_bytes = 0;
 
  std::size_t m_size = 0;
 
  bool m_sorted = false;
 
  std::array<void*, s_numColumns> m_colPtrs{};
 
  std::vector<std::size_t> m_cellStart;
 
  Box m_sortBox;
 
  RealVect m_sortDx = RealVect::Zero;
 
  RealVect m_sortProbLo = RealVect::Zero;
 
  std::size_t
  cellKey(const std::size_t a_index, const Box& a_box, const RealVect& a_dx, const RealVect& a_probLo) const noexcept
  {
    const IntVect  loEnd   = a_box.smallEnd();
    const IntVect  boxSize = a_box.size();
    const RealVect x       = this->position(a_index);
 
    std::size_t lin    = 0;
    std::size_t stride = 1;
    for (int d = 0; d < SpaceDim; d++) {
      int c = static_cast<int>(std::floor((x[d] - a_probLo[d]) / a_dx[d])) - loEnd[d];
      if (c < 0) {
        c = 0;
      }
      if (c >= boxSize[d]) {
        c = boxSize[d] - 1;
      }
      lin += static_cast<std::size_t>(c) * stride;
      stride *= static_cast<std::size_t>(boxSize[d]);
    }
 
    return lin;
  }
 
  bool
  checkCellSort(const Box& a_box, const RealVect& a_dx, const RealVect& a_probLo) const noexcept
  {
    if (!m_sorted) {
      return false;
    }
    if (m_cellStart.size() != static_cast<std::size_t>(a_box.numPts()) + 1) {
      return false;
    }
    if (m_cellStart.back() != m_size) {
      return false;
    }
 
    for (std::size_t c = 0; c < this->numCells(); c++) {
      if (m_cellStart[c] > m_cellStart[c + 1]) {
        return false;
      }
 
      for (std::size_t i = m_cellStart[c]; i < m_cellStart[c + 1]; i++) {
        if (this->cellKey(i, a_box, a_dx, a_probLo) != c) {
          return false;
        }
      }
    }
 
    return true;
  }
 
  static std::size_t
  alignUp(const std::size_t a_x) noexcept
  {
    return (a_x + (s_align - 1)) / s_align * s_align;
  }
 
  void
  steal(ParticleSoA& a_o) noexcept
  {
    m_buffer       = a_o.m_buffer;
    m_capacity     = a_o.m_capacity;
    m_bytes        = a_o.m_bytes;
    m_size         = a_o.m_size;
    m_sorted       = a_o.m_sorted;
    m_sortBox      = a_o.m_sortBox;
    m_sortDx       = a_o.m_sortDx;
    m_sortProbLo   = a_o.m_sortProbLo;
    m_colPtrs      = a_o.m_colPtrs;
    m_cellStart    = std::move(a_o.m_cellStart);
    a_o.m_buffer   = nullptr;
    a_o.m_capacity = 0;
    a_o.m_bytes    = 0;
    a_o.m_size     = 0;
    a_o.m_sorted   = false;
    a_o.m_colPtrs  = {};
  }
 
  void
  growTo(const std::size_t a_needed)
  {
    if (a_needed <= m_capacity) {
      return;
    }
 
    // Saturating geometric growth: 2 * m_capacity overflows std::size_t once m_capacity exceeds
    // SIZE_MAX/2, and reserve()'s "a_capacity <= m_capacity" guard would then silently no-op.
    const std::size_t doubled = (m_capacity == 0) ? 16
                                                  : (m_capacity <= std::numeric_limits<std::size_t>::max() / 2
                                                       ? 2 * m_capacity
                                                       : std::numeric_limits<std::size_t>::max());
 
    this->reserve(std::max(a_needed, doubled));
  }
 
  void
  growIfFull()
  {
    if (m_size == m_capacity) {
      this->growTo(m_size + 1);
    }
  }
 
  template <std::size_t... I>
  std::size_t
  computeOffsets(const std::size_t                      a_capacity,
                 std::array<std::size_t, s_numColumns>& a_off,
                 std::index_sequence<I...>) const noexcept
  {
    std::size_t running = 0;
    using expander      = int[];
    (void)expander{0, (a_off[I] = running, running += alignUp(a_capacity * sizeof(ColumnType<I>)), 0)...};
    return running;
  }
 
  template <std::size_t... I>
  void
  setColPtrs(const std::array<std::size_t, s_numColumns>& a_off, std::index_sequence<I...>) noexcept
  {
    using expander = int[];
    (void)expander{0, (m_colPtrs[I] = static_cast<char*>(m_buffer) + a_off[I], 0)...};
  }
 
  template <std::size_t... I>
  void
  moveColumns(void* a_newBuf, const std::array<std::size_t, s_numColumns>& a_off, std::index_sequence<I...>) noexcept
  {
    using expander = int[];
    (void)expander{
      0,
      (std::memcpy(static_cast<char*>(a_newBuf) + a_off[I], m_colPtrs[I], m_size * sizeof(ColumnType<I>)), 0)...};
  }
 
  template <std::size_t... I>
  void
  copyLiveColumnsTo(ParticleSoA& a_dst, std::index_sequence<I...>) const noexcept
  {
    using expander = int[];
    (void)expander{0, (std::memcpy(a_dst.m_colPtrs[I], m_colPtrs[I], m_size * sizeof(ColumnType<I>)), 0)...};
  }
 
  template <std::size_t... I>
  void
  appendColumnsFrom(const ParticleSoA& a_other, const std::size_t a_at, std::index_sequence<I...>) noexcept
  {
    using expander = int[];
    (void)expander{0,
                   (std::memcpy(static_cast<char*>(m_colPtrs[I]) + a_at * sizeof(ColumnType<I>),
                                a_other.m_colPtrs[I],
                                a_other.m_size * sizeof(ColumnType<I>)),
                    0)...};
  }
 
  template <std::size_t... I>
  void
  copyParticleColumns(const ParticleSoA& a_src,
                      const std::size_t  a_srcIdx,
                      const std::size_t  a_dstIdx,
                      std::index_sequence<I...>) noexcept
  {
    using expander = int[];
    (void)expander{0,
                   (this->template columnByIndex<I>()[a_dstIdx] = a_src.template columnByIndex<I>()[a_srcIdx], 0)...};
  }
 
  void
  reallocate(const std::size_t a_capacity);
 
  void
  permute(const std::vector<std::size_t>& a_dest);
 
  template <std::size_t K>
  void
  scatterColumn(void* a_newBuf, const std::size_t a_off, const std::vector<std::size_t>& a_dest) noexcept
  {
    const ColumnType<K>* src = this->template columnByIndex<K>();
    ColumnType<K>*       dst = reinterpret_cast<ColumnType<K>*>(static_cast<char*>(a_newBuf) + a_off);
    for (std::size_t i = 0; i < m_size; i++) {
      dst[a_dest[i]] = src[i];
    }
  }
 
  template <std::size_t... I>
  void
  scatterColumns(void*                                        a_newBuf,
                 const std::array<std::size_t, s_numColumns>& a_off,
                 const std::vector<std::size_t>&              a_dest,
                 std::index_sequence<I...>) noexcept
  {
    using expander = int[];
    (void)expander{0, ((void)this->template scatterColumn<I>(a_newBuf, a_off[I], a_dest), 0)...};
  }
 
  template <std::size_t... J>
  void
  writePayload(const P& a_p, std::index_sequence<J...>) noexcept
  {
    constexpr auto cols = Traits::columns;
    using expander      = int[];
    (void)expander{0,
                   ((void)(this->template columnByIndex<s_payloadBegin + J>()[m_size] = a_p.*std::get<J>(cols)), 0)...};
  }
 
  template <std::size_t... J>
  P
  gatherPayload(const std::size_t a_index, std::index_sequence<J...>) const
  {
    constexpr auto cols = Traits::columns;
    P              p;
    using expander = int[];
    (void)expander{0,
                   ((void)(p.*std::get<J>(cols) = this->template columnByIndex<s_payloadBegin + J>()[a_index]), 0)...};
    return p;
  }
 
  template <std::size_t... I>
  void
  removeImpl(const std::size_t a_index, std::index_sequence<I...>) noexcept
  {
    using expander = int[];
    (void)expander{
      0,
      ((void)(this->template columnByIndex<I>()[a_index] = this->template columnByIndex<I>()[m_size - 1]), 0)...};
  }
 
  template <std::size_t... I>
  void
  linearizeImpl(unsigned char* a_buf, const std::size_t a_index, std::index_sequence<I...>) const noexcept
  {
    using expander = int[];
    (void)expander{0, ((void)detail::pushBytes(a_buf, this->template columnByIndex<I>()[a_index]), 0)...};
  }
 
  template <std::size_t... I>
  void
  delinearizeImpl(const unsigned char* a_buf, std::index_sequence<I...>) noexcept
  {
    using expander = int[];
    (void)expander{0, ((void)this->template pullColumn<I>(a_buf), 0)...};
  }
 
  template <std::size_t K>
  void
  pullColumn(const unsigned char*& a_buf) noexcept
  {
    detail::pullBytes(a_buf, this->template columnByIndex<K>()[m_size]);
  }
 
  template <std::size_t... J>
  static constexpr std::size_t
  h5PayloadSize(std::index_sequence<J...>) noexcept
  {
    return detail::sumSizes(sizeof(ColumnType<s_payloadBegin + J>)...);
  }
 
  template <std::size_t... J>
  void
  h5LinearizePayload(unsigned char*& a_buf, const std::size_t a_index, std::index_sequence<J...>) const noexcept
  {
    using expander = int[];
    (void)expander{0,
                   ((void)detail::pushBytes(a_buf, this->template columnByIndex<s_payloadBegin + J>()[a_index]), 0)...};
  }
 
  template <std::size_t... J>
  void
  h5DelinearizePayload(const unsigned char*& a_buf, std::index_sequence<J...>) noexcept
  {
    using expander = int[];
    (void)expander{0,
                   ((void)detail::pullBytes(a_buf, this->template columnByIndex<s_payloadBegin + J>()[m_size]), 0)...};
  }
 
  template <std::size_t... I>
  std::size_t
  byteSpanImpl(std::index_sequence<I...>) const noexcept
  {
    std::size_t span = 0;
    using expander   = int[];
    (void)expander{0,
                   (span = (span > (static_cast<std::size_t>(static_cast<const char*>(m_colPtrs[I]) -
                                                             static_cast<const char*>(m_buffer)) +
                                    m_size * sizeof(ColumnType<I>)))
                             ? span
                             : (static_cast<std::size_t>(static_cast<const char*>(m_colPtrs[I]) -
                                                         static_cast<const char*>(m_buffer)) +
                                m_size * sizeof(ColumnType<I>)),
                    0)...};
    return span;
  }
};
 
#include <CD_NamespaceFooter.H>
 
#include <CD_ParticleSoAImplem.H>
 
#endif