smash::Grid< Options > Class Template Reference

#include <grid.h>

template<GridOptions Options = GridOptions::Normal>
class smash::Grid< Options >

Abstracts a list of cells that partition the particles in the experiment into regions of space that can interact / cannot interact.

This class is used to construct a helper data structure to reduce the combinatorics of finding particle pairs that could interact (scatter). It takes a list of ParticleData objects and sorts them in such a way that it is easy to look only at lists of particles that have a chance of interacting.

Template Parameters

Options

This policy parameter determines whether ghost cells are created to support periodic boundaries, or not.

Definition at line 79 of file grid.h.

Inheritance diagram for smash::Grid< Options >:

Collaboration diagram for smash::Grid< Options >:

Public Member Functions
	Grid (const Particles &particles, double min_cell_length, double timestep_duration, CellSizeStrategy strategy=CellSizeStrategy::Optimal)
	Constructs a grid from the given particle list particles. More...
	Grid (const std::pair< std::array< double, 3 >, std::array< double, 3 >> &min_and_length, const Particles &particles, double min_cell_length, double timestep_duration, CellSizeStrategy strategy=CellSizeStrategy::Optimal)
	Constructs a grid with the given minimum grid coordinates and grid length. More...
void	iterate_cells (const std::function< void(const ParticleList &)> &search_cell_callback, const std::function< void(const ParticleList &, const ParticleList &)> &neighbor_cell_callback) const
	Iterates over all cells in the grid and calls the callback arguments with a search cell and 0 to 13 neighbor cells. More...
double	cell_volume () const
template<>
void	iterate_cells (const std::function< void(const ParticleList &)> &search_cell_callback, const std::function< void(const ParticleList &, const ParticleList &)> &neighbor_cell_callback) const
	Specialization of iterate_cells. More...
template<>
void	iterate_cells (const std::function< void(const ParticleList &)> &search_cell_callback, const std::function< void(const ParticleList &, const ParticleList &)> &neighbor_cell_callback) const
	Specialization of iterate_cells. More...

Private Member Functions
SizeType	make_index (SizeType x, SizeType y, SizeType z) const
SizeType	make_index (std::array< SizeType, 3 > idx) const

Private Attributes
const std::array< double, 3 >	length_
	The 3 lengths of the complete grid. Used for periodic boundary wrapping. More...
double	cell_volume_
	The volume of a single cell. More...
std::array< int, 3 >	number_of_cells_
	The number of cells in x, y, and z direction. More...
std::vector< ParticleList >	cells_
	The cell storage. More...

Additional Inherited Members
Public Types inherited from smash::GridBase
typedef int	SizeType
	A type to store the sizes. More...
Static Protected Member Functions inherited from smash::GridBase
static std::pair< std::array< double, 3 >, std::array< double, 3 > >	find_min_and_length (const Particles &particles)

Constructor & Destructor Documentation

Grid() [1/2]

template<GridOptions Options = GridOptions::Normal>

smash::Grid< Options >::Grid ( const Particles &  particles, double  min_cell_length, double  timestep_duration, CellSizeStrategy  strategy = CellSizeStrategy::Optimal  )

inline

Constructs a grid from the given particle list particles.

It automatically determines the necessary size for the grid from the positions of the particles.

Parameters

[in]	particles	The particles to place onto the grid.
[in]	min_cell_length	The minimal length a cell must have.
[in]	timestep_duration	Duration of the timestep. It is necessary for formation times treatment: if particle is fully or partially formed before the end of the timestep, it has to be on the grid.
[in]	strategy	The strategy for determining the cell size

Definition at line 93 of file grid.h.

      : Grid{find_min_and_length(particles), std::move(particles),
             min_cell_length, timestep_duration, strategy} {}

Grid() [2/2]

template<GridOptions O>

template smash::Grid< Options >::Grid	(	const std::pair< std::array< double, 3 >, std::array< double, 3 >> &	min_and_length,
		const Particles &	particles,
		double	min_cell_length,
		double	timestep_duration,
		CellSizeStrategy	strategy = CellSizeStrategy::Optimal
	)

Constructs a grid with the given minimum grid coordinates and grid length.

If you need periodic boundaries you have to use this constructor to set the correct length to use for wrapping particles around the borders.

Parameters

[in]	min_and_length	A pair consisting of the three min coordinates and the three lengths.
[in]	particles	The particles to place onto the grid.
[in]	min_cell_length	The minimal length a cell must have.
[in]	timestep_duration	duration of the timestep in fm/c
[in]	strategy	The strategy for determining the cell size

Exceptions

runtime_error

if your box length is smaller than the grid length.

Definition at line 101 of file grid.cc.

    : length_(min_and_length.second) {
  const auto min_position = min_and_length.first;
  const SizeType particle_count = particles.size();
 
  // very simple setup for non-periodic boundaries and largest cellsize strategy
  if (O == GridOptions::Normal && strategy == CellSizeStrategy::Largest) {
    number_of_cells_ = {1, 1, 1};
    cell_volume_ = length_[0] * length_[1] * length_[2];
    cells_.clear();
    cells_.reserve(1);
    cells_.emplace_back(particles.copy_to_vector());
    return;
  }
 
  // The number of cells is determined by the min and max coordinates where
  // particles are positioned and the maximal interaction length (which equals
  // the length of a cell).
  // But don't let the number of cells exceed the actual number of particles.
  // That would be overkill. Let max_cells³ ≤ particle_count (conversion to
  // int truncates).
  // Consider that particle placement into cells uses half-open intervals. Thus
  // a cell includes particles in [0, a[. The next cell [a, 2a[. And so on. This
  // is important for calculating the number of cells. If length * index_factor
  // (equivalent to length / max_interaction_length) is integral, then
  // length * index_factor + 1 determines the number of required cells. That's
  // because the last cell will then store particles in the interval
  // [length, length + max_interaction_length[. The code below achieves this
  // effect by rounding down (floor) and adding 1 afterwards.
  const int max_cells =
      (O == GridOptions::Normal)
          ? std::cbrt(particle_count)
          : std::max(2, static_cast<int>(std::cbrt(particle_count)));
 
  std::string error_box_too_small =
      "Input error: Your box is too small for the grid.\n"
      "The minimal length of the box is given by: " +
      std::to_string(2 * max_interaction_length) +
      " fm with the given timestep size.\n"
      "If you have large timesteps please reduce them.\n"
      "A larger box or the use of testparticles also helps.\n"
      "Please take a look at your config.";
 
  // This normally equals 1/max_interaction_length, but if the number of cells
  // is reduced (because of low density) then this value is smaller.
  std::array<double, 3> index_factor = {1. / max_interaction_length,
                                        1. / max_interaction_length,
                                        1. / max_interaction_length};
  for (std::size_t i = 0; i < number_of_cells_.size(); ++i) {
    number_of_cells_[i] =
        (strategy == CellSizeStrategy::Largest)
            ? 2
            : static_cast<int>(std::floor(length_[i] * index_factor[i])) +
                  // The last cell in each direction can be smaller than
                  // max_interaction_length. In that case periodic boundaries
                  // will not work correctly. Thus, we need to reduce the number
                  // of cells in that direction by one and make the last cell
                  // larger. This basically merges a smaller boundary cell into
                  // a full cell inside the grid. There's a ~0% chance that the
                  // given boundaries create an integral number of cells with
                  // length of max_interaction_length. Therefore, just make the
                  // default number of cells one less than for non-periodic
                  // boundaries.
                  (O == GridOptions::Normal ? 1 : 0);
    // Only in the case of periodic boundaries (i.e. GridOptions != Normal) the
    // number of cells can be zero.
    if (number_of_cells_[i] == 0) {
      // The minimal cell length exceeds the length of the box.
      throw std::runtime_error(error_box_too_small);
    }
    // std::nextafter implements a safety margin so that no valid position
    // inside the grid can reference an out-of-bounds cell
    if (number_of_cells_[i] > max_cells) {
      number_of_cells_[i] = max_cells;
      index_factor[i] = number_of_cells_[i] / length_[i];
      while (index_factor[i] * length_[i] >= number_of_cells_[i]) {
        index_factor[i] = std::nextafter(index_factor[i], 0.);
      }
      assert(index_factor[i] * length_[i] < number_of_cells_[i]);
    } else if (O == GridOptions::PeriodicBoundaries) {
      if (number_of_cells_[i] == 1) {
        number_of_cells_[i] = 2;
      }
      index_factor[i] = number_of_cells_[i] / length_[i];
      while (index_factor[i] * length_[i] >= number_of_cells_[i]) {
        index_factor[i] = std::nextafter(index_factor[i], 0.);
      }
      assert(index_factor[i] * length_[i] < number_of_cells_[i]);
      // Verify that cell length did not become smaller than
      // the max. interaction length by increasing the number of cells
      if (1. / index_factor[i] <= std::nextafter(max_interaction_length, 0.)) {
        // The minimal cell length exceeds
        // the length of a grid cell in the box.
        throw std::runtime_error(error_box_too_small);
      }
    }
  }
 
  cell_volume_ = (length_[0] / number_of_cells_[0]) *
                 (length_[1] / number_of_cells_[1]) *
                 (length_[2] / number_of_cells_[2]);
 
  if (O == GridOptions::Normal &&
      all_of(number_of_cells_, [](SizeType n) { return n <= 2; })) {
    // dilute limit:
    // the grid would have <= 2x2x2 cells, meaning every particle has to be
    // compared with every other particle anyway. Then we can just as well
    // fall back to not using the grid at all
    // For a grid with periodic boundaries the situation is different and we
    // never want to have a grid smaller than 2x2x2.
    logg[LGrid].debug(
        "There would only be ", number_of_cells_,
        " cells. Therefore the Grid falls back to a single cell / "
        "particle list.");
    number_of_cells_ = {1, 1, 1};
    cell_volume_ = length_[0] * length_[1] * length_[2];
    cells_.resize(1);
    cells_.front().reserve(particles.size());
    std::copy_if(particles.begin(), particles.end(),
                 std::back_inserter(cells_.front()),
                 [&](const ParticleData &p) {
                   return p.xsec_scaling_factor(timestep_duration) > 0.0;
                 });  // filter out the particles that can not interact
  } else {
    // construct a normal grid
    logg[LGrid].debug("min: ", min_position, "\nlength: ", length_,
                      "\ncell_volume: ", cell_volume_,
                      "\ncells: ", number_of_cells_,
                      "\nindex_factor: ", index_factor);
 
    // After the grid parameters are determined, we can start placing the
    // particles in cells.
    cells_.resize(number_of_cells_[0] * number_of_cells_[1] *
                  number_of_cells_[2]);
 
    // Returns the one-dimensional cell-index from the position vector inside
    // the grid.
    // This simply calculates the distance to min_position and multiplies it
    // with index_factor to determine the 3 x,y,z indexes to pass to make_index.
    auto &&cell_index_for = [&](const ParticleData &p) {
      return make_index(
          std::floor((p.position()[1] - min_position[0]) * index_factor[0]),
          std::floor((p.position()[2] - min_position[1]) * index_factor[1]),
          std::floor((p.position()[3] - min_position[2]) * index_factor[2]));
    };
 
    for (const auto &p : particles) {
      if (p.xsec_scaling_factor(timestep_duration) > 0.0) {
        const auto idx = cell_index_for(p);
#ifndef NDEBUG
        if (idx >= SizeType(cells_.size())) {
          logg[LGrid].fatal(
              source_location,
              "\nan out-of-bounds access would be necessary for the "
              "particle ",
              p, "\nfor a grid with the following parameters:\nmin: ",
              min_position, "\nlength: ", length_,
              "\ncells: ", number_of_cells_, "\nindex_factor: ", index_factor,
              "\ncells_.size: ", cells_.size(), "\nrequested index: ", idx);
          throw std::runtime_error("out-of-bounds grid access on construction");
        }
#endif
        cells_[idx].push_back(p);
      }
    }
  }
 
  logg[LGrid].debug(cells_);
}

Member Function Documentation

iterate_cells() [1/3]

template<GridOptions Options = GridOptions::Normal>

void smash::Grid< Options >::iterate_cells	(	const std::function< void(const ParticleList &)> &	search_cell_callback,
		const std::function< void(const ParticleList &, const ParticleList &)> &	neighbor_cell_callback
	)		const

Iterates over all cells in the grid and calls the callback arguments with a search cell and 0 to 13 neighbor cells.

The neighbor cells are constructed like this:

• one cell at x+1
• three cells (x-1,x,x+1) at y+1
• nine cells (x-1, y-1)...(x+1, y+1) at z+1

Parameters

[in]	search_cell_callback	A callable called for/with every non-empty cell in the grid.
[in]	neighbor_cell_callback	A callable called for/with every non-empty cell and adjacent cell combination. For a periodic grid, the first argument will be adjusted to wrap around the grid.

cell_volume()

template<GridOptions Options = GridOptions::Normal>

double smash::Grid< Options >::cell_volume ( ) const

inline

Returns

the volume of a single grid cell

Definition at line 142 of file grid.h.

{ return cell_volume_; }

make_index() [1/2]

template<GridOptions Options>

Grid< Options >::SizeType smash::Grid< Options >::make_index ( SizeType  x, SizeType  y, SizeType  z  ) const

inlineprivate

Returns

the one-dimensional cell-index from the 3-dim index x, y, z.

Definition at line 275 of file grid.cc.

                                              {
  return (z * number_of_cells_[1] + y) * number_of_cells_[0] + x;
}

Here is the caller graph for this function:

make_index() [2/2]

template<GridOptions Options = GridOptions::Normal>

SizeType smash::Grid< Options >::make_index ( std::array< SizeType, 3 >  idx ) const

inlineprivate

Returns

the one-dimensional cell-index from the 3-dim index idx. This is a convenience overload for the above function.

Definition at line 155 of file grid.h.

                                                       {
    return make_index(idx[0], idx[1], idx[2]);
  }

iterate_cells() [2/3]

template<>

void smash::Grid< GridOptions::Normal >::iterate_cells	(	const std::function< void(const ParticleList &)> &	search_cell_callback,
		const std::function< void(const ParticleList &, const ParticleList &)> &	neighbor_cell_callback
	)		const

Specialization of iterate_cells.

Definition at line 288 of file grid.cc.

                                       {
  std::array<SizeType, 3> search_index;
  SizeType &x = search_index[0];
  SizeType &y = search_index[1];
  SizeType &z = search_index[2];
  SizeType search_cell_index = 0;
  for (z = 0; z < number_of_cells_[2]; ++z) {
    for (y = 0; y < number_of_cells_[1]; ++y) {
      for (x = 0; x < number_of_cells_[0]; ++x, ++search_cell_index) {
        assert(search_cell_index == make_index(search_index));
        assert(search_cell_index >= 0);
        assert(search_cell_index < SizeType(cells_.size()));
        const ParticleList &search = cells_[search_cell_index];
        search_cell_callback(search);
 
        const auto &dz_list = z == number_of_cells_[2] - 1 ? ZERO : ZERO_ONE;
        const auto &dy_list = number_of_cells_[1] == 1
                                  ? ZERO
                                  : y == 0 ? ZERO_ONE
                                           : y == number_of_cells_[1] - 1
                                                 ? MINUS_ONE_ZERO
                                                 : MINUS_ONE_ZERO_ONE;
        const auto &dx_list = number_of_cells_[0] == 1
                                  ? ZERO
                                  : x == 0 ? ZERO_ONE
                                           : x == number_of_cells_[0] - 1
                                                 ? MINUS_ONE_ZERO
                                                 : MINUS_ONE_ZERO_ONE;
        for (SizeType dz : dz_list) {
          for (SizeType dy : dy_list) {
            for (SizeType dx : dx_list) {
              const auto di = make_index(dx, dy, dz);
              if (di > 0) {
                neighbor_cell_callback(search, cells_[search_cell_index + di]);
              }
            }
          }
        }
      }
    }
  }
}

iterate_cells() [3/3]

template<>

void smash::Grid< GridOptions::PeriodicBoundaries >::iterate_cells	(	const std::function< void(const ParticleList &)> &	search_cell_callback,
		const std::function< void(const ParticleList &, const ParticleList &)> &	neighbor_cell_callback
	)		const

Specialization of iterate_cells.

Definition at line 354 of file grid.cc.

                                       {
  std::array<SizeType, 3> search_index;
  SizeType &x = search_index[0];
  SizeType &y = search_index[1];
  SizeType &z = search_index[2];
  SizeType search_cell_index = 0;
 
  // defaults:
  std::array<NeighborLookup, 2> dz_list;
  std::array<NeighborLookup, 3> dy_list;
  std::array<NeighborLookup, 3> dx_list;
 
  assert(number_of_cells_[2] >= 2);
  assert(number_of_cells_[1] >= 2);
  assert(number_of_cells_[0] >= 2);
 
  for (z = 0; z < number_of_cells_[2]; ++z) {
    dz_list[0].index = z;
    dz_list[1].index = z + 1;
    if (dz_list[1].index == number_of_cells_[2]) {
      dz_list[1].index = 0;
      // last z in the loop, so no need to reset wrap again
      dz_list[1].wrap = NeedsToWrap::MinusLength;
    }
    for (y = 0; y < number_of_cells_[1]; ++y) {
      dy_list[0].index = y;
      dy_list[1].index = y - 1;
      dy_list[2].index = y + 1;
      dy_list[2].wrap = NeedsToWrap::No;
      if (y == 0) {
        dy_list[1] = dy_list[2];
        dy_list[2].index = number_of_cells_[1] - 1;
        dy_list[2].wrap = NeedsToWrap::PlusLength;
      } else if (dy_list[2].index == number_of_cells_[1]) {
        dy_list[2].index = 0;
        dy_list[2].wrap = NeedsToWrap::MinusLength;
      }
      for (x = 0; x < number_of_cells_[0]; ++x, ++search_cell_index) {
        dx_list[0].index = x;
        dx_list[1].index = x - 1;
        dx_list[2].index = x + 1;
        dx_list[2].wrap = NeedsToWrap::No;
        if (x == 0) {
          dx_list[1] = dx_list[2];
          dx_list[2].index = number_of_cells_[0] - 1;
          dx_list[2].wrap = NeedsToWrap::PlusLength;
        } else if (dx_list[2].index == number_of_cells_[0]) {
          dx_list[2].index = 0;
          dx_list[2].wrap = NeedsToWrap::MinusLength;
        }
 
        assert(search_cell_index == make_index(search_index));
        assert(search_cell_index >= 0);
        assert(search_cell_index < SizeType(cells_.size()));
        ParticleList search = cells_[search_cell_index];
        search_cell_callback(search);
 
        auto virtual_search_index = search_index;
        ThreeVector wrap_vector = {};  // no change
        auto current_wrap_vector = wrap_vector;
 
        for (const auto &dz : dz_list) {
          if (dz.wrap == NeedsToWrap::MinusLength) {
            // last dz in the loop, so no need to undo the wrap
            wrap_vector[2] = -length_[2];
            virtual_search_index[2] = -1;
          }
          for (const auto &dy : dy_list) {
            // only the last dy in dy_list can wrap
            if (dy.wrap == NeedsToWrap::MinusLength) {
              wrap_vector[1] = -length_[1];
              virtual_search_index[1] = -1;
            } else if (dy.wrap == NeedsToWrap::PlusLength) {
              wrap_vector[1] = length_[1];
              virtual_search_index[1] = number_of_cells_[1];
            }
            for (const auto &dx : dx_list) {
              // only the last dx in dx_list can wrap
              if (dx.wrap == NeedsToWrap::MinusLength) {
                wrap_vector[0] = -length_[0];
                virtual_search_index[0] = -1;
              } else if (dx.wrap == NeedsToWrap::PlusLength) {
                wrap_vector[0] = length_[0];
                virtual_search_index[0] = number_of_cells_[0];
              }
              assert(dx.index >= 0);
              assert(dx.index < number_of_cells_[0]);
              assert(dy.index >= 0);
              assert(dy.index < number_of_cells_[1]);
              assert(dz.index >= 0);
              assert(dz.index < number_of_cells_[2]);
              const auto neighbor_cell_index =
                  make_index(dx.index, dy.index, dz.index);
              assert(neighbor_cell_index >= 0);
              assert(neighbor_cell_index < SizeType(cells_.size()));
              if (neighbor_cell_index <= make_index(virtual_search_index)) {
                continue;
              }
 
              if (wrap_vector != current_wrap_vector) {
                logg[LGrid].debug("translating search cell by ",
                                  wrap_vector - current_wrap_vector);
                for_each(search, [&](ParticleData &p) {
                  p = p.translated(wrap_vector - current_wrap_vector);
                });
                current_wrap_vector = wrap_vector;
              }
              neighbor_cell_callback(search, cells_[neighbor_cell_index]);
            }
            virtual_search_index[0] = search_index[0];
            wrap_vector[0] = 0;
          }
          virtual_search_index[1] = search_index[1];
          wrap_vector[1] = 0;
        }
      }
    }
  }
}

Member Data Documentation

length_

template<GridOptions Options = GridOptions::Normal>

const std::array<double, 3> smash::Grid< Options >::length_

private

The 3 lengths of the complete grid. Used for periodic boundary wrapping.

Definition at line 160 of file grid.h.

cell_volume_

template<GridOptions Options = GridOptions::Normal>

double smash::Grid< Options >::cell_volume_

private

The volume of a single cell.

Definition at line 163 of file grid.h.

number_of_cells_

template<GridOptions Options = GridOptions::Normal>

std::array<int, 3> smash::Grid< Options >::number_of_cells_

private

The number of cells in x, y, and z direction.

Definition at line 166 of file grid.h.

cells_

template<GridOptions Options = GridOptions::Normal>

std::vector<ParticleList> smash::Grid< Options >::cells_

private

The cell storage.

Definition at line 169 of file grid.h.

---

The documentation for this class was generated from the following files:

• src/include/smash/grid.h
• src/grid.cc