/******************************************************************************* * This file is part of SWIFT. * Copyright (c) 2024 Will Roper (w.roper@sussex.ac.uk) * * This program is free software: you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published * by the Free Software Foundation, either version 3 of the License, or * (at your option) any later version. * * This program is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU General Public License for more details. * * You should have received a copy of the GNU Lesser General Public License * along with this program. If not, see . * ******************************************************************************/ /* Config parameters. */ #include /* Local headers. */ #include "active.h" #include "cell.h" #include "engine.h" #include "gravity_properties.h" #include "runner.h" #include "runner_doiact_grav.h" #include "space.h" #include "timers.h" /** * @brief Performs M-M interactions between a given top-level cell and * all other top level cells not interacted with via pair tasks. * * This is the non-periodic case where there is no mesh so all cells not * handled by a pair task are interacted with here in this long range * gravity function. * * This function is used when running a non-periodic uniform box (i.e. non-zoom * simulation) and just loops over every other cell with particles and interacts * with them. * * @param r The thread #runner. * @param ci The #cell of interest. * @param top The top-level parent of the #cell of interest. */ void runner_do_grav_long_range_uniform_non_periodic(struct runner *r, struct cell *ci, struct cell *top) { struct engine *e = r->e; struct space *s = e->s; /* Get the mutlipole of the cell we are interacting. */ struct gravity_tensors *const multi_i = ci->grav.multipole; /* Recover the list of top-level cells */ struct cell *cells = e->s->cells_top; int *cells_with_particles = e->s->cells_with_particles_top; const int nr_cells_with_particles = e->s->nr_cells_with_particles; /* Loop over all the top-level cells and go for a M-M interaction if * well-separated */ for (int n = 0; n < nr_cells_with_particles; ++n) { /* Handle on the top-level cell and it's gravity business*/ struct cell *cj = &cells[cells_with_particles[n]]; struct gravity_tensors *const multi_j = cj->grav.multipole; /* Avoid self contributions */ if (top == cj) continue; /* Skip empty cells */ if (multi_j->m_pole.M_000 == 0.f) continue; if (cell_can_use_pair_mm(top, cj, e, e->s, /*use_rebuild_data=*/1, /*is_tree_walk=*/0, /*periodic boundaries*/ s->periodic, /*use_mesh*/ s->periodic)) { /* Call the PM interaction function on the active sub-cells of ci */ runner_dopair_grav_mm_nonsym(r, ci, cj); // runner_dopair_recursive_grav_pm(r, ci, cj); /* Record that this multipole received a contribution */ multi_i->pot.interacted = 1; } /* We are in charge of this pair */ } /* Loop over top-level cells */ } /** * @brief Performs M-M interactions between a given top-level cell and * all other top level cells not interacted with via pair tasks. * * This is the non-periodic case where there is no mesh so all cells not * handled by a pair task are interacted with here in this long range * gravity function. * * This function is used when running a non-periodic zoom simulation and * will loop over all non-zoom cells but use the void cell hierarchy to * interact with all zoom cells. * * @param r The thread #runner. * @param ci The #cell of interest. * @param top The top-level parent of the #cell of interest. */ void runner_do_grav_long_range_zoom_non_periodic(struct runner *r, struct cell *ci, struct cell *top) { struct engine *e = r->e; struct space *s = e->s; /* Get the mutlipole of the cell we are interacting. */ struct gravity_tensors *const multi_i = ci->grav.multipole; /* Recover the list of top-level cells */ struct cell *bkg_cells = e->s->zoom_props->bkg_cells_top; #ifdef SWIFT_DEBUG_CHECKS /* Define counters used to count gparts. */ size_t tested_gparts = 0; #endif /* Since the zoom cells will be handled by the void cell hierarchy we can * just loop over all other cells which are not zoom cells. This is * trivial since the zoom cells are first in cells_top. */ for (int cjd = 0; cjd < s->zoom_props->nr_bkg_cells; cjd++) { /* Handle on the top-level cell and it's gravity business*/ struct cell *cj = &bkg_cells[cjd]; struct gravity_tensors *const multi_j = cj->grav.multipole; /* Skip empty cells */ if (multi_j->m_pole.M_000 == 0.f) continue; #ifdef SWIFT_DEBUG_CHECKS tested_gparts += multi_j->m_pole.num_gpart; #endif /* Avoid self contributions */ if (top == cj) continue; if (cell_can_use_pair_mm(top, cj, e, e->s, /*use_rebuild_data=*/1, /*is_tree_walk=*/0, /*periodic boundaries*/ s->periodic, /*use_mesh*/ s->periodic)) { /* Call the PM interaction function on the active sub-cells of ci */ runner_dopair_grav_mm_nonsym(r, ci, cj); // runner_dopair_recursive_grav_pm(r, ci, cj); /* Record that this multipole received a contribution */ multi_i->pot.interacted = 1; } /* We are in charge of this pair */ } /* Loop over top-level cells */ #ifdef SWIFT_DEBUG_CHECKS /* Ensure we at leasted against all possible gparts. */ if (tested_gparts != e->s->nr_gparts) { error( "Not all gparts were tested in long range gravity task! (tested: %ld, " "total: %ld)", tested_gparts, e->s->nr_gparts); } #endif } /** * @brief Performs M-M interactions between a given top-level cell and all other * top level cells not interacted with via pair tasks or the mesh. * * This is used when running a zoom simulation, the space is periodic and there * is a mesh, therefore we only interact with cells that are closer than the * mesh interaction distance but further than the direct interaction distance. * * This function will be "clever" and only loop over the parts of the * top-level grid that are not covered by the mesh. Add some buffer for * safety. * * This function handles the following long range interactions: * - zoom -> bkg * - zoom -> buffer (if their are buffer cells) * - bkg -> bkg * - bkg -> buffer (if their are buffer cells) * * Since we define the original pair tasks at the void level we only need to * consider combinations of background cells. However, If a zoom cell has a long * range task it means there were no tasks in the void level but we still only * need to consider interactions between the zoom cell and the background cells. * Note that the top cell pointer here is already going to be the top level void * parent cell if ci is a zoom cell. * * @param r The thread #runner. * @param ci The #cell of interest. * @param top The top-level parent of the #cell of interest. */ void runner_do_grav_long_range_zoom_periodic(struct runner *r, struct cell *ci, struct cell *top) { struct engine *e = r->e; struct space *s = e->s; struct cell *bkg_cells = s->zoom_props->bkg_cells_top; const int periodic = e->mesh->periodic; const double dim[3] = {e->mesh->dim[0], e->mesh->dim[1], e->mesh->dim[2]}; /* Get the maximum distance at which we can have a non-mesh interaction. */ const double max_distance = e->mesh->r_cut_max; const double max_distance2 = max_distance * max_distance; /* Get the mutlipole of the cell we are interacting. */ struct gravity_tensors *const multi_i = ci->grav.multipole; /* Get the (i,j,k) location of the top-level cell in the grid. */ int top_i = top->loc[0] * s->iwidth[0]; int top_j = top->loc[1] * s->iwidth[1]; int top_k = top->loc[2] * s->iwidth[2]; /* Maximal distance any interaction can take place before the mesh kicks in, * rounded up to the next integer */ int d = ceil(max_distance * max3(s->iwidth[0], s->iwidth[1], s->iwidth[2])) + 1; /* Loop over plausibly useful cells */ for (int ii = top_i - d; ii <= top_i + d; ++ii) { for (int jj = top_j - d; jj <= top_j + d; ++jj) { for (int kk = top_k - d; kk <= top_k + d; ++kk) { /* Box wrap */ const int iii = (ii + s->cdim[0]) % s->cdim[0]; const int jjj = (jj + s->cdim[1]) % s->cdim[1]; const int kkk = (kk + s->cdim[2]) % s->cdim[2]; /* Get the cell */ const int cell_index = cell_getid(s->cdim, iii, jjj, kkk); /* Handle on the top-level cell */ struct cell *cj = &bkg_cells[cell_index]; /* Avoid self contributions */ if (top == cj) continue; /* Handle on the top-level cell's gravity business*/ const struct gravity_tensors *multi_j = cj->grav.multipole; /* Skip empty cells */ if (multi_j->m_pole.M_000 == 0.f) continue; /* Minimal distance between any pair of particles */ const double min_radius2 = cell_min_dist2(top, cj, periodic, dim); /* Are we beyond the distance where the truncated forces are 0 ?*/ if (min_radius2 > max_distance2) { /* Record that this multipole received a contribution */ multi_i->pot.interacted = 1; /* We are done here. */ continue; } /* Shall we interact with this cell? */ if (cell_can_use_pair_mm(top, cj, e, e->s, /*use_rebuild_data=*/1, /*is_tree_walk=*/0, /*periodic boundaries*/ s->periodic, /*use_mesh*/ s->periodic)) { /* Call the PM interaction function on the active sub-cells of ci */ runner_dopair_grav_mm_nonsym(r, ci, cj); /* Record that this multipole received a contribution */ multi_i->pot.interacted = 1; } /* We can interact with this cell */ } /* Loop over relevant top-level cells (k) */ } /* Loop over relevant top-level cells (j) */ } /* Loop over relevant top-level cells (i) */ } /** * @brief Performs M-M interactions between a given top-level cell and all other * top level cells not interacted with via pair tasks or the mesh. * * This is used when running a uniform box, the space is periodic and there is a * mesh, therefore we only interact with cells that are closer than the mesh * interaction distance but further than the direct interaction distance. * * This function will be "clever" and only loop over the parts of the * top-level grid that are not covered by the mesh. Add some buffer for * safety. * * @param r The thread #runner. * @param ci The #cell of interest. * @param top The top-level parent of the #cell of interest. */ void runner_do_grav_long_range_uniform_periodic(struct runner *r, struct cell *ci, struct cell *top) { struct engine *e = r->e; struct space *s = e->s; struct cell *cells = s->cells_top; const int periodic = e->mesh->periodic; const double dim[3] = {e->mesh->dim[0], e->mesh->dim[1], e->mesh->dim[2]}; /* Get the maximum distance at which we can have a non-mesh interaction. */ const double max_distance = e->mesh->r_cut_max; const double max_distance2 = max_distance * max_distance; /* Get the mutlipole of the cell we are interacting. */ struct gravity_tensors *const multi_i = ci->grav.multipole; /* Get the (i,j,k) location of the top-level cell in the grid. */ int top_i = top->loc[0] * s->iwidth[0]; int top_j = top->loc[1] * s->iwidth[1]; int top_k = top->loc[2] * s->iwidth[2]; /* Maximal distance any interaction can take place before the mesh kicks in, * rounded up to the next integer */ int d = ceil(max_distance * max3(s->iwidth[0], s->iwidth[1], s->iwidth[2])) + 1; /* Loop over plausibly useful cells */ for (int ii = top_i - d; ii <= top_i + d; ++ii) { for (int jj = top_j - d; jj <= top_j + d; ++jj) { for (int kk = top_k - d; kk <= top_k + d; ++kk) { /* Box wrap */ const int iii = (ii + s->cdim[0]) % s->cdim[0]; const int jjj = (jj + s->cdim[1]) % s->cdim[1]; const int kkk = (kk + s->cdim[2]) % s->cdim[2]; /* Get the cell */ const int cell_index = cell_getid(s->cdim, iii, jjj, kkk); /* Handle on the top-level cell */ struct cell *cj = &cells[cell_index]; /* Avoid self contributions */ if (top == cj) continue; /* Handle on the top-level cell's gravity business*/ const struct gravity_tensors *multi_j = cj->grav.multipole; /* Skip empty cells */ if (multi_j->m_pole.M_000 == 0.f) continue; /* Minimal distance between any pair of particles */ const double min_radius2 = cell_min_dist2(top, cj, periodic, dim); /* Are we beyond the distance where the truncated forces are 0 ?*/ if (min_radius2 > max_distance2) { /* Record that this multipole received a contribution */ multi_i->pot.interacted = 1; /* We are done here. */ continue; } /* Shall we interact with this cell? */ if (cell_can_use_pair_mm(top, cj, e, e->s, /*use_rebuild_data=*/1, /*is_tree_walk=*/0, /*periodic boundaries*/ s->periodic, /*use_mesh*/ s->periodic)) { /* Call the PM interaction function on the active sub-cells of ci */ runner_dopair_grav_mm_nonsym(r, ci, cj); // runner_dopair_recursive_grav_pm(r, ci, cj); /* Record that this multipole received a contribution */ multi_i->pot.interacted = 1; } /* We can interact with this cell */ } /* Loop over relevant top-level cells (k) */ } /* Loop over relevant top-level cells (j) */ } /* Loop over relevant top-level cells (i) */ } /** * @brief Accumalate the number of particle mesh interactions for debugging * purposes. * * This is the varaint used when running a zoom simulation. * * @param r The thread #runner. * @param ci The #cell of interest. * @param top The current top-level cell. */ void runner_count_mesh_interactions_zoom(struct runner *r, struct cell *ci, struct cell *top) { #if defined(SWIFT_DEBUG_CHECKS) || defined(SWIFT_GRAVITY_FORCE_CHECKS) struct engine *e = r->e; struct space *s = e->s; struct cell *cells = s->cells_top; const int periodic = e->mesh->periodic; const double dim[3] = {e->mesh->dim[0], e->mesh->dim[1], e->mesh->dim[2]}; /* Get the maximum distance at which we can have a non-mesh interaction. */ const double max_distance = e->mesh->r_cut_max; const double max_distance2 = max_distance * max_distance; /* Get the multipole of the cell we are interacting. */ struct gravity_tensors *const multi_i = ci->grav.multipole; /* Loop over all cells. */ for (int n = 0; n < s->nr_cells; n++) { /* Skip void cells to avoid double counting their top level progeny * in the zoom (and buffer) cell grids. */ if (cells[n].subtype == cell_subtype_void) continue; /* Handle on the top-level cell and it's gravity business*/ struct cell *cj = &cells[n]; struct gravity_tensors *const multi_j = cj->grav.multipole; /* Get the top level cell of the current cj */ struct cell *top_j = cj->top; /* If we are in a zoom cell we need to jump up the void hierarchy * to get the top level cell. */ if (top_j->void_parent != NULL) { top_j = cj->void_parent->top; } /* If we had buffer cells then we may need an extra jump since the * top-level for a zoom cell is at: * zoom->top->void_parent->top->void_parent->top. */ if (top_j->void_parent != NULL) { top_j = top_j->void_parent->top; } /* Avoid self contributions */ if (top == top_j) continue; /* Skip empty cells */ if (multi_j->m_pole.M_000 == 0.f) continue; /* Minimal distance between any pair of particles */ const double min_radius2 = cell_min_dist2(top, top_j, periodic, dim); /* Are we beyond the distance where the truncated forces are 0 ?*/ if (min_radius2 > max_distance2) { #ifdef SWIFT_DEBUG_CHECKS /* Need to account for the interactions we missed */ accumulate_add_ll(&multi_i->pot.num_interacted, multi_j->m_pole.num_gpart); #endif #ifdef SWIFT_GRAVITY_FORCE_CHECKS /* Need to account for the interactions we missed */ accumulate_add_ll(&multi_i->pot.num_interacted_pm, multi_j->m_pole.num_gpart); #endif /* Record that this multipole received a contribution */ multi_i->pot.interacted = 1; } } #else error( "This function should not be called without debugging checks or " "force checks enabled!"); #endif } /** * @brief Accumalate the number of particle mesh interactions for debugging * purposes. * * This is the varaint used when running a uniform box simulation. * * @param r The thread #runner. * @param ci The #cell of interest. * @param top The current top-level cell. */ void runner_count_mesh_interactions_uniform(struct runner *r, struct cell *ci, struct cell *top) { #if defined(SWIFT_DEBUG_CHECKS) || defined(SWIFT_GRAVITY_FORCE_CHECKS) struct engine *e = r->e; struct space *s = e->s; struct cell *cells = s->cells_top; const int periodic = e->mesh->periodic; const double dim[3] = {e->mesh->dim[0], e->mesh->dim[1], e->mesh->dim[2]}; /* Get the maximum distance at which we can have a non-mesh interaction. */ const double max_distance = e->mesh->r_cut_max; const double max_distance2 = max_distance * max_distance; /* Get the multipole of the cell we are interacting. */ struct gravity_tensors *const multi_i = ci->grav.multipole; /* Loop over all cells. */ for (int n = 0; n < s->nr_cells; n++) { /* Handle on the top-level cell and it's gravity business*/ struct cell *cj = &cells[n]; struct gravity_tensors *const multi_j = cj->grav.multipole; /* Avoid self contributions */ if (top == cj) continue; /* Skip empty cells */ if (multi_j->m_pole.M_000 == 0.f) continue; /* Minimal distance between any pair of particles */ const double min_radius2 = cell_min_dist2(top, cj, periodic, dim); /* Are we beyond the distance where the truncated forces are 0 ?*/ if (min_radius2 > max_distance2) { #ifdef SWIFT_DEBUG_CHECKS /* Need to account for the interactions we missed */ accumulate_add_ll(&multi_i->pot.num_interacted, multi_j->m_pole.num_gpart); #endif #ifdef SWIFT_GRAVITY_FORCE_CHECKS /* Need to account for the interactions we missed */ accumulate_add_ll(&multi_i->pot.num_interacted_pm, multi_j->m_pole.num_gpart); #endif /* Record that this multipole received a contribution */ multi_i->pot.interacted = 1; } } #else error( "This function should not be called without debugging checks or " "force checks enabled!"); #endif } /** * @brief Performs all M-M interactions between a given top-level cell and * all the other top-levels that are far enough. * * This function is the main long-range gravity function. It is called by * the runner and will call the appropriate interaction function for the * given cell type/space periodicity. * * @param r The thread #runner. * @param ci The #cell of interest. * @param timer Are we timing this ? */ void runner_do_grav_long_range(struct runner *r, struct cell *ci, const int timer) { TIMER_TIC; struct space *s = r->e->s; /* Is the space periodic? */ const int periodic = s->periodic; /* Anything to do here? */ if (!cell_is_active_gravity(ci, r->e)) return; if (ci->nodeID != engine_rank) error("Non-local cell in long-range gravity task!"); /* Check multipole has been drifted */ if (ci->grav.ti_old_multipole < r->e->ti_current) cell_drift_multipole(ci, r->e); /* Find this cell's top-level (great-)parent */ struct cell *top = ci; while (top->parent != NULL) top = top->parent; /* If we have a nested cell the true top level cell where we defined the * interactions is the background void parent cell. */ /* TODO: This is a bit of a hack, we should actually have the top pointers set * properly during void tree splitting. */ while (top->void_parent != NULL) top = top->void_parent->top; /* Call the appropriate interaction function based on the type of the * cell in question. */ if (periodic) { switch (ci->type) { case cell_type_regular: runner_do_grav_long_range_uniform_periodic(r, ci, top); break; case cell_type_zoom: runner_do_grav_long_range_zoom_periodic(r, ci, top); break; case cell_type_buffer: runner_do_grav_long_range_zoom_periodic(r, ci, top); break; case cell_type_bkg: runner_do_grav_long_range_zoom_periodic(r, ci, top); break; default: error("Unknown cell type in long-range gravity task!"); } } else { switch (ci->type) { case cell_type_regular: runner_do_grav_long_range_uniform_non_periodic(r, ci, top); break; case cell_type_zoom: runner_do_grav_long_range_zoom_non_periodic(r, ci, top); break; case cell_type_buffer: runner_do_grav_long_range_zoom_non_periodic(r, ci, top); break; case cell_type_bkg: runner_do_grav_long_range_zoom_non_periodic(r, ci, top); break; default: error("Unknown cell type in long-range gravity task!"); } } #if defined(SWIFT_DEBUG_CHECKS) || defined(SWIFT_GRAVITY_FORCE_CHECKS) /* Count the number of mesh interactions if using the mesh. */ if (periodic && s->with_zoom_region) { runner_count_mesh_interactions_zoom(r, ci, top); } else if (periodic) { runner_count_mesh_interactions_uniform(r, ci, top); } #endif if (timer) TIMER_TOC(timer_dograv_long_range); }