22 static std::ostream &operator<<(std::ostream &out,
const std::vector<T> &v) {
23 auto column = out.tellp();
25 for (
const auto &x : v) {
26 if (out.tellp() - column >= 100) {
36 static std::ostream &operator<<(std::ostream &out,
37 const std::initializer_list<T> &v) {
38 auto column = out.tellp();
40 for (
const auto &x : v) {
41 if (out.tellp() - column >= 100) {
50 template <
typename T, std::
size_t N>
51 static std::ostream &operator<<(std::ostream &out,
const std::array<T, N> &a) {
52 auto column = out.tellp();
54 for (
const auto &x : a) {
55 if (out.tellp() - column >= 100) {
66 static constexpr
int LGrid = LogArea::Grid::id;
71 std::pair<std::array<double, 3>, std::array<double, 3>>
73 std::pair<std::array<double, 3>, std::array<double, 3>> r;
74 auto &min_position = r.first;
75 auto &length = r.second;
80 min_position = {{first_position[1], first_position[2], first_position[3]}};
81 auto max_position = min_position;
82 for (
const auto &
p : particles) {
83 const auto &pos =
p.position();
84 min_position[0] = std::min(min_position[0], pos[1]);
85 min_position[1] = std::min(min_position[1], pos[2]);
86 min_position[2] = std::min(min_position[2], pos[3]);
87 max_position[0] = std::max(max_position[0], pos[1]);
88 max_position[1] = std::max(max_position[1], pos[2]);
89 max_position[2] = std::max(max_position[2], pos[3]);
91 length[0] = max_position[0] - min_position[0];
92 length[1] = max_position[1] - min_position[1];
93 length[2] = max_position[2] - min_position[2];
100 template <Gr
idOptions O>
103 const Particles &particles,
double max_interaction_length,
105 : length_(min_and_length.second) {
106 const auto min_position = min_and_length.first;
110 if (O == GridOptions::Normal && strategy == CellSizeStrategy::Largest) {
111 number_of_cells_ = {1, 1, 1};
112 cell_volume_ = length_[0] * length_[1] * length_[2];
133 const int max_cells =
134 (O == GridOptions::Normal)
135 ? std::cbrt(particle_count)
136 : std::max(2, static_cast<int>(std::cbrt(particle_count)));
138 std::string error_box_too_small =
139 "Input error: Your box is too small for the grid.\n"
140 "The minimal length of the box is given by: " +
141 std::to_string(2 * max_interaction_length) +
142 " fm with the given timestep size.\n"
143 "If you have large timesteps please reduce them.\n"
144 "A larger box or the use of testparticles also helps.\n"
145 "Please take a look at your config.";
149 std::array<double, 3> index_factor = {1. / max_interaction_length,
150 1. / max_interaction_length,
151 1. / max_interaction_length};
152 for (std::size_t i = 0; i < number_of_cells_.size(); ++i) {
153 number_of_cells_[i] =
154 (strategy == CellSizeStrategy::Largest)
156 : static_cast<int>(std::floor(length_[i] * index_factor[i])) +
167 (O == GridOptions::Normal ? 1 : 0);
170 if (number_of_cells_[i] == 0) {
172 throw std::runtime_error(error_box_too_small);
176 if (number_of_cells_[i] > max_cells) {
177 number_of_cells_[i] = max_cells;
178 index_factor[i] = number_of_cells_[i] / length_[i];
179 while (index_factor[i] * length_[i] >= number_of_cells_[i]) {
180 index_factor[i] = std::nextafter(index_factor[i], 0.);
182 assert(index_factor[i] * length_[i] < number_of_cells_[i]);
183 }
else if (O == GridOptions::PeriodicBoundaries) {
184 if (number_of_cells_[i] == 1) {
185 number_of_cells_[i] = 2;
187 index_factor[i] = number_of_cells_[i] / length_[i];
188 while (index_factor[i] * length_[i] >= number_of_cells_[i]) {
189 index_factor[i] = std::nextafter(index_factor[i], 0.);
191 assert(index_factor[i] * length_[i] < number_of_cells_[i]);
194 if (1. / index_factor[i] <= std::nextafter(max_interaction_length, 0.)) {
197 throw std::runtime_error(error_box_too_small);
202 cell_volume_ = (length_[0] / number_of_cells_[0]) *
203 (length_[1] / number_of_cells_[1]) *
204 (length_[2] / number_of_cells_[2]);
206 if (O == GridOptions::Normal &&
215 "There would only be ", number_of_cells_,
216 " cells. Therefore the Grid falls back to a single cell / "
218 number_of_cells_ = {1, 1, 1};
219 cell_volume_ = length_[0] * length_[1] * length_[2];
221 cells_.front().reserve(particles.
size());
222 std::copy_if(particles.
begin(), particles.
end(),
223 std::back_inserter(cells_.front()),
225 return p.xsec_scaling_factor(timestep_duration) > 0.0;
229 logg[
LGrid].debug(
"min: ", min_position,
"\nlength: ", length_,
230 "\ncell_volume: ", cell_volume_,
231 "\ncells: ", number_of_cells_,
232 "\nindex_factor: ", index_factor);
236 cells_.resize(number_of_cells_[0] * number_of_cells_[1] *
237 number_of_cells_[2]);
245 std::floor((
p.position()[1] - min_position[0]) * index_factor[0]),
246 std::floor((
p.position()[2] - min_position[1]) * index_factor[1]),
247 std::floor((
p.position()[3] - min_position[2]) * index_factor[2]));
250 for (
const auto &
p : particles) {
251 if (
p.xsec_scaling_factor(timestep_duration) > 0.0) {
252 const auto idx = cell_index_for(
p);
254 if (idx >=
SizeType(cells_.size())) {
257 "\nan out-of-bounds access would be necessary for the "
259 p,
"\nfor a grid with the following parameters:\nmin: ",
260 min_position,
"\nlength: ", length_,
261 "\ncells: ", number_of_cells_,
"\nindex_factor: ", index_factor,
262 "\ncells_.size: ", cells_.size(),
"\nrequested index: ", idx);
263 throw std::runtime_error(
"out-of-bounds grid access on construction");
266 cells_[idx].push_back(
p);
274 template <Gr
idOptions Options>
277 return (z * number_of_cells_[1] + y) * number_of_cells_[0] + x;
280 static const std::initializer_list<GridBase::SizeType>
ZERO{0};
281 static const std::initializer_list<GridBase::SizeType>
ZERO_ONE{0, 1};
289 const std::function<
void(
const ParticleList &)> &search_cell_callback,
290 const std::function<
void(
const ParticleList &,
const ParticleList &)>
291 &neighbor_cell_callback)
const {
292 std::array<SizeType, 3> search_index;
297 for (z = 0; z < number_of_cells_[2]; ++z) {
298 for (y = 0; y < number_of_cells_[1]; ++y) {
299 for (x = 0; x < number_of_cells_[0]; ++x, ++search_cell_index) {
300 assert(search_cell_index == make_index(search_index));
301 assert(search_cell_index >= 0);
302 assert(search_cell_index <
SizeType(cells_.size()));
303 const ParticleList &search = cells_[search_cell_index];
304 search_cell_callback(search);
306 const auto &dz_list = z == number_of_cells_[2] - 1 ?
ZERO :
ZERO_ONE;
307 const auto &dy_list = number_of_cells_[1] == 1
310 : y == number_of_cells_[1] - 1
313 const auto &dx_list = number_of_cells_[0] == 1
316 : x == number_of_cells_[0] - 1
322 const auto di = make_index(dx, dy, dz);
324 neighbor_cell_callback(search, cells_[search_cell_index + di]);
355 const std::function<
void(
const ParticleList &)> &search_cell_callback,
356 const std::function<
void(
const ParticleList &,
const ParticleList &)>
357 &neighbor_cell_callback)
const {
358 std::array<SizeType, 3> search_index;
365 std::array<NeighborLookup, 2> dz_list;
366 std::array<NeighborLookup, 3> dy_list;
367 std::array<NeighborLookup, 3> dx_list;
369 assert(number_of_cells_[2] >= 2);
370 assert(number_of_cells_[1] >= 2);
371 assert(number_of_cells_[0] >= 2);
373 for (z = 0; z < number_of_cells_[2]; ++z) {
374 dz_list[0].index = z;
375 dz_list[1].index = z + 1;
376 if (dz_list[1].index == number_of_cells_[2]) {
377 dz_list[1].index = 0;
379 dz_list[1].wrap = NeedsToWrap::MinusLength;
381 for (y = 0; y < number_of_cells_[1]; ++y) {
382 dy_list[0].index = y;
383 dy_list[1].index = y - 1;
384 dy_list[2].index = y + 1;
385 dy_list[2].wrap = NeedsToWrap::No;
387 dy_list[1] = dy_list[2];
388 dy_list[2].index = number_of_cells_[1] - 1;
389 dy_list[2].wrap = NeedsToWrap::PlusLength;
390 }
else if (dy_list[2].index == number_of_cells_[1]) {
391 dy_list[2].index = 0;
392 dy_list[2].wrap = NeedsToWrap::MinusLength;
394 for (x = 0; x < number_of_cells_[0]; ++x, ++search_cell_index) {
395 dx_list[0].index = x;
396 dx_list[1].index = x - 1;
397 dx_list[2].index = x + 1;
398 dx_list[2].wrap = NeedsToWrap::No;
400 dx_list[1] = dx_list[2];
401 dx_list[2].index = number_of_cells_[0] - 1;
402 dx_list[2].wrap = NeedsToWrap::PlusLength;
403 }
else if (dx_list[2].index == number_of_cells_[0]) {
404 dx_list[2].index = 0;
405 dx_list[2].wrap = NeedsToWrap::MinusLength;
408 assert(search_cell_index == make_index(search_index));
409 assert(search_cell_index >= 0);
410 assert(search_cell_index <
SizeType(cells_.size()));
411 ParticleList search = cells_[search_cell_index];
412 search_cell_callback(search);
414 auto virtual_search_index = search_index;
416 auto current_wrap_vector = wrap_vector;
418 for (
const auto &dz : dz_list) {
419 if (dz.wrap == NeedsToWrap::MinusLength) {
421 wrap_vector[2] = -length_[2];
422 virtual_search_index[2] = -1;
424 for (
const auto &dy : dy_list) {
426 if (dy.wrap == NeedsToWrap::MinusLength) {
427 wrap_vector[1] = -length_[1];
428 virtual_search_index[1] = -1;
429 }
else if (dy.wrap == NeedsToWrap::PlusLength) {
430 wrap_vector[1] = length_[1];
431 virtual_search_index[1] = number_of_cells_[1];
433 for (
const auto &dx : dx_list) {
435 if (dx.wrap == NeedsToWrap::MinusLength) {
436 wrap_vector[0] = -length_[0];
437 virtual_search_index[0] = -1;
438 }
else if (dx.wrap == NeedsToWrap::PlusLength) {
439 wrap_vector[0] = length_[0];
440 virtual_search_index[0] = number_of_cells_[0];
442 assert(dx.index >= 0);
443 assert(dx.index < number_of_cells_[0]);
444 assert(dy.index >= 0);
445 assert(dy.index < number_of_cells_[1]);
446 assert(dz.index >= 0);
447 assert(dz.index < number_of_cells_[2]);
448 const auto neighbor_cell_index =
449 make_index(dx.index, dy.index, dz.index);
450 assert(neighbor_cell_index >= 0);
451 assert(neighbor_cell_index <
SizeType(cells_.size()));
452 if (neighbor_cell_index <= make_index(virtual_search_index)) {
456 if (wrap_vector != current_wrap_vector) {
457 logg[
LGrid].debug(
"translating search cell by ",
458 wrap_vector - current_wrap_vector);
460 p =
p.translated(wrap_vector - current_wrap_vector);
462 current_wrap_vector = wrap_vector;
464 neighbor_cell_callback(search, cells_[neighbor_cell_index]);
466 virtual_search_index[0] = search_index[0];
469 virtual_search_index[1] = search_index[1];
478 const std::pair<std::array<double, 3>, std::array<double, 3>>
480 const Particles &particles,
double max_interaction_length,
483 const std::pair<std::array<double, 3>, std::array<double, 3>>
485 const Particles &particles,
double max_interaction_length,