/******************************************************************************* * This file is part of SWIFT. * Copyright (c) 2019 Stuart McAlpine (stuart.mcalpine@helsinki.fi) * 2024 Will J. 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 */ #include /* Includes */ #include /* Local includes */ #include "cell.h" #include "engine.h" #include "gravity_properties.h" #include "space.h" #include "zoom.h" /* mpi headers. */ #ifdef WITH_MPI #include #endif /* Declare the task diff grav constant. */ int zoom_bkg_subdepth_diff_grav = zoom_bkg_subdepth_diff_grav_default; /** * @brief Read parameter file for "ZoomRegion" properties. * * @param params Swift parameter structure. * @param props The zoom properties struct. */ void zoom_parse_params(struct swift_params *params, struct zoom_region_properties *props) { /* Set the zoom cdim. */ props->zoom_cell_depth = parser_get_opt_param_int(params, "ZoomRegion:zoom_top_level_depth", 2); /* Set the target background cdim, default is a negative value so that if no * value is given for a target then the zoom region defines the background * cell size. */ int bkg_cdim = parser_get_opt_param_int(params, "ZoomRegion:bkg_top_level_cells", space_max_top_level_cells_default); for (int i = 0; i < 3; i++) { props->bkg_cdim[i] = bkg_cdim; } /* Get the ratio between the zoom region size and buffer cell size. * Ignored if buffer cells aren't needed. */ props->buffer_cell_depth = parser_get_opt_param_int(params, "ZoomRegion:buffer_top_level_depth", 0); /* Ensure the buffer cell depth is less than the zoom cell depth. */ if (props->buffer_cell_depth > props->zoom_cell_depth) { error("Buffer cell depth must be less than the zoom cell depth."); } /* Extract the zoom width boost factor (used to define the buffer around the * zoom region). */ props->region_pad_factor = parser_get_opt_param_float(params, "ZoomRegion:region_pad_factor", 1.1); /* Extract the depth we'll split neighbour cells to. */ props->neighbour_max_tree_depth = parser_get_opt_param_int( params, "ZoomRegion:neighbour_max_tree_depth", -1); /* Extract the minimum difference between the task level and the leaves * for background cells. */ zoom_bkg_subdepth_diff_grav = parser_get_opt_param_int(params, "ZoomRegion:bkg_subdepth_diff_grav", zoom_bkg_subdepth_diff_grav_default); } /** * @brief Compute the zoom region centre and boundaries. * * Finds the dimensions of the high resolution particle distribution and * computes the necessary shift to shift the zoom region to the centre of the * box. This shift is stored to be applied in space_init and for * transformation when writing out. * * @param s The space */ double zoom_get_region_dim_and_shift(struct space *s) { /* Initialise values we will need. */ const size_t nr_gparts = s->nr_gparts; double min_bounds[3] = {FLT_MAX, FLT_MAX, FLT_MAX}; double max_bounds[3] = {-FLT_MAX, -FLT_MAX, -FLT_MAX}; double midpoint[3] = {0.0, 0.0, 0.0}; double com[3] = {0.0, 0.0, 0.0}; double mtot = 0.0; double ini_dims[3] = {0.0, 0.0, 0.0}; const double box_mid[3] = {s->dim[0] / 2.0, s->dim[1] / 2.0, s->dim[2] / 2.0}; /* Find the min/max location in each dimension for each * high resolution gravity particle (non-background), and their COM. */ for (size_t k = 0; k < nr_gparts; k++) { /* Skip background particles. */ if (s->gparts[k].type == swift_type_dark_matter_background) { continue; } /* Unpack the particle positions. * NOTE: these will have already been shifted by the user requested amount * in space_init if shift in the parameter file is non-zero. */ const double x = s->gparts[k].x[0]; const double y = s->gparts[k].x[1]; const double z = s->gparts[k].x[2]; /* Wrap if periodic. */ if (s->periodic) { box_wrap(x, 0.0, s->dim[0]); box_wrap(y, 0.0, s->dim[1]); box_wrap(z, 0.0, s->dim[2]); } /* Ammend boundaries for this particle. */ if (x > max_bounds[0]) max_bounds[0] = x; if (y > max_bounds[1]) max_bounds[1] = y; if (z > max_bounds[2]) max_bounds[2] = z; if (x < min_bounds[0]) min_bounds[0] = x; if (y < min_bounds[1]) min_bounds[1] = y; if (z < min_bounds[2]) min_bounds[2] = z; /* Total up mass and position for COM. */ mtot += s->gparts[k].mass; com[0] += x * s->gparts[k].mass; com[1] += y * s->gparts[k].mass; com[2] += z * s->gparts[k].mass; } #ifdef WITH_MPI /* Share answers amoungst nodes. */ /* Boundary. */ MPI_Allreduce(MPI_IN_PLACE, &min_bounds[0], 3, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD); MPI_Allreduce(MPI_IN_PLACE, &max_bounds[1], 3, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD); /* CoM. */ MPI_Allreduce(MPI_IN_PLACE, com, 3, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); MPI_Allreduce(MPI_IN_PLACE, &mtot, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); #endif /* Finalize CoM calcuation. */ const double imass = 1.0 / mtot; com[0] *= imass; com[1] *= imass; com[2] *= imass; /* Get the initial dimensions and midpoint. */ for (int i = 0; i < 3; i++) { ini_dims[i] = max_bounds[i] - min_bounds[i]; midpoint[i] = min_bounds[i] + (ini_dims[i] / 2.0); } /* Calculate the shift needed to place the mid point of the high res * particles at the centre of the box. This shift is applied to the * particles in space_init in space.c */ for (int i = 0; i < 3; i++) { s->zoom_props->zoom_shift[i] = box_mid[i] - midpoint[i]; } /* We shouldn't shift if the shift is incremental. */ for (int i = 0; i < 3; i++) { if (fabs(s->zoom_props->zoom_shift[i]) < 0.01 * s->dim[i]) { s->zoom_props->zoom_shift[i] = 0.0; } } /* If the volume isn't periodic then we can't shift. */ if (!s->periodic) { for (int i = 0; i < 3; i++) { if (fabs(s->zoom_props->zoom_shift[i]) > 0.01 * s->dim[i]) { error( "Cannot shift the zoom region to the centre of the box " "when the box is not periodic. Centre the CoM of the high " "resolution particles in the box. (shift=[%f, %f, %f], dim=[%f, " "%f, %f])", s->zoom_props->zoom_shift[0], s->zoom_props->zoom_shift[1], s->zoom_props->zoom_shift[2], s->dim[0], s->dim[1], s->dim[2]); } } } /* Let's shift the COM. * NOTE: boundaries are recalculated relative to box centre later. */ for (int i = 0; i < 3; i++) s->zoom_props->com[i] += s->zoom_props->zoom_shift[i]; /* Compute maximum side length of the zoom region, we need zoom dim to be * equal. */ double ini_dim = max3(ini_dims[0], ini_dims[1], ini_dims[2]); return ini_dim; } /** * @brief Compute the void region geometry. * * The void region is the region covered by background cells above the zoom * region. If the void region is sufficiently close to the zoom region size * then the two will be made equivalent later on. Otherwide the void region is * equivalent to the buffer region. * * This function will derive the bounds of the void region and return how many * zoom regions tesselate the void region. * * @param s The space * @param ini_max_dim The dim of the zoom region before tesselating the * volume. * * @return The number of zoom regions that tesselate the void region. */ int zoom_get_void_geometry(struct space *s, const double region_dim) { /* Get the lower and upper bounds of the zoom region based on the initial * dimensions and the padding factor. */ double lower_bounds[3]; double upper_bounds[3]; for (int i = 0; i < 3; i++) { lower_bounds[i] = (s->dim[i] / 2) - (region_dim / 2.0); upper_bounds[i] = (s->dim[i] / 2) + (region_dim / 2.0); } /* Assign these temporary bounds to the zoom region bounds, we'll overwrite * these later but useful to have them for now. */ for (int i = 0; i < 3; i++) { s->zoom_props->region_lower_bounds[i] = lower_bounds[i]; s->zoom_props->region_upper_bounds[i] = upper_bounds[i]; } /* Find the background cell edges that contain these bounds. */ double void_lower_bounds[3]; double void_upper_bounds[3]; for (int i = 0; i < 3; i++) { int lower = (int)floor(lower_bounds[i] * s->iwidth[i]); int upper = (int)floor(upper_bounds[i] * s->iwidth[i]); void_lower_bounds[i] = lower * s->width[i]; void_upper_bounds[i] = (upper + 1) * s->width[i]; } /* Assign the void bounds. */ for (int i = 0; i < 3; i++) { s->zoom_props->void_lower_bounds[i] = void_lower_bounds[i]; s->zoom_props->void_upper_bounds[i] = void_upper_bounds[i]; } /* Compute the void region dimensions. */ for (int i = 0; i < 3; i++) { s->zoom_props->void_dim[i] = void_upper_bounds[i] - void_lower_bounds[i]; } /* Compute the number of zoom regions that tesselate the void region. */ int nr_zoom_regions = ceil(s->zoom_props->void_dim[0] / region_dim); return nr_zoom_regions; } /** * @brief Compute the number of child cells in a single parent cell at a given * depth. * * @param region_dim The dimension of the region. * @param parent_width The width of the parent cell. * @param child_depth The depth of the child cell within the parent. * * @return The number of child cells in a single parent cell. */ static int zoom_get_cdim_at_depth(double region_dim, double parent_width, int child_depth) { /* How many parent_widths are in the region? (ensure correct rounding) */ int region_parent_cdim = floor((region_dim + (0.1 * parent_width)) / parent_width); /* We now know how many parent cells we have in the region, use this and the * depth of the zoom region to calculate the cdim (the number of parents times * the number of children in a parent. */ return region_parent_cdim * pow(2, child_depth); } void zoom_get_geometry_no_buffer_cells(struct space *s) { /* If we have a buffer cell depth warn that we will ignore it. */ if (s->zoom_props->buffer_cell_depth > 0) { message("No buffer cells are needed, ignoring buffer cell depth."); s->zoom_props->buffer_cell_depth = 0; } /* Zero the buffer region properties explictly. */ for (int i = 0; i < 3; i++) { s->zoom_props->buffer_lower_bounds[i] = 0.0; s->zoom_props->buffer_upper_bounds[i] = 0.0; s->zoom_props->buffer_dim[i] = 0.0; s->zoom_props->buffer_cdim[i] = 0; s->zoom_props->buffer_width[i] = 0.0; } /* Match the zoom reigon bounds to the void region bounds. */ for (int i = 0; i < 3; i++) { s->zoom_props->region_lower_bounds[i] = s->zoom_props->void_lower_bounds[i]; s->zoom_props->region_upper_bounds[i] = s->zoom_props->void_upper_bounds[i]; } /* Compute the zoom region dimensions. */ for (int i = 0; i < 3; i++) { s->zoom_props->dim[i] = s->zoom_props->region_upper_bounds[i] - s->zoom_props->region_lower_bounds[i]; } /* Compute the number of zoom cells in the void region. */ int cdim = zoom_get_cdim_at_depth(s->zoom_props->dim[0], s->width[0], s->zoom_props->zoom_cell_depth); /* Compute the zoom cdim and cell width. */ for (int i = 0; i < 3; i++) { s->zoom_props->cdim[i] = cdim; s->zoom_props->width[i] = s->zoom_props->dim[i] / cdim; s->zoom_props->iwidth[i] = 1.0 / s->zoom_props->width[i]; } } /** * @brief Compute the geometry of the zoom region with buffer cells. * * This function computes the geometry of the zoom region when buffer cells are * enabled. It calculates the bounds, dimensions, and cell widths for both the * buffer and zoom regions. * * Currently, buffer cells are not fully supported and this function will * simply through an error if called. * * @param s The space */ void zoom_get_geometry_with_buffer_cells(struct space *s) { error( "Buffer cells currently provide no performance benefit and carry a " "significant complexity cost. They are thus not currently fully " "supported. Set ZoomRegion:buffer_top_level_depth to 0 to disable " "buffer cells."); /* Ensure we have a buffer cell depth. */ if (s->zoom_props->buffer_cell_depth == 0) { error( "Current cell structure requires buffer cells but not buffer cell has " "been given. ZoomRegion:buffer_top_level_depth must be greater than " "0."); } /* Match the buffer region bounds to the void region bounds. */ for (int i = 0; i < 3; i++) { s->zoom_props->buffer_lower_bounds[i] = s->zoom_props->void_lower_bounds[i]; s->zoom_props->buffer_upper_bounds[i] = s->zoom_props->void_upper_bounds[i]; } /* Compute the buffer region dimensions. */ for (int i = 0; i < 3; i++) { s->zoom_props->buffer_dim[i] = s->zoom_props->buffer_upper_bounds[i] - s->zoom_props->buffer_lower_bounds[i]; } /* Compute the number of buffer cells in the void region. */ int buffer_cdim = zoom_get_cdim_at_depth(s->zoom_props->buffer_dim[0], s->width[0], s->zoom_props->buffer_cell_depth); /* Compute the buffer cdim and cell width. */ for (int i = 0; i < 3; i++) { s->zoom_props->buffer_cdim[i] = buffer_cdim; s->zoom_props->buffer_width[i] = s->zoom_props->buffer_dim[i] / buffer_cdim; s->zoom_props->buffer_iwidth[i] = 1.0 / s->zoom_props->buffer_width[i]; } /* Find the buffer cell edges that contain the zoom region bounds. */ double region_lower_bounds[3]; double region_upper_bounds[3]; for (int i = 0; i < 3; i++) { int lower = (int)floor((s->zoom_props->region_lower_bounds[i] - s->zoom_props->buffer_lower_bounds[i]) * s->zoom_props->buffer_iwidth[i]); int upper = (int)floor((s->zoom_props->region_upper_bounds[i] - s->zoom_props->buffer_lower_bounds[i]) * s->zoom_props->buffer_iwidth[i]); region_lower_bounds[i] = lower * s->zoom_props->buffer_width[i] + s->zoom_props->buffer_lower_bounds[i]; region_upper_bounds[i] = (upper + 1) * s->zoom_props->buffer_width[i] + s->zoom_props->buffer_lower_bounds[i]; } /* Assign the new aligned zoom bounds. */ for (int i = 0; i < 3; i++) { s->zoom_props->region_lower_bounds[i] = region_lower_bounds[i]; s->zoom_props->region_upper_bounds[i] = region_upper_bounds[i]; } /* Compute the zoom region dimensions. */ for (int i = 0; i < 3; i++) { s->zoom_props->dim[i] = region_upper_bounds[i] - region_lower_bounds[i]; } /* Compute the number of zoom cells in the zoom region. Here we need to * subtract the buffer depth from the user defined zoom depth, both are * defined from the background cells but the calculation below is from the * buffer cells down to the zoom level. */ int cdim = zoom_get_cdim_at_depth( s->zoom_props->dim[0], s->zoom_props->buffer_width[0], s->zoom_props->zoom_cell_depth - s->zoom_props->buffer_cell_depth); /* Compute the zoom cdim and cell width. */ for (int i = 0; i < 3; i++) { s->zoom_props->cdim[i] = cdim; s->zoom_props->width[i] = s->zoom_props->dim[i] / cdim; s->zoom_props->iwidth[i] = 1.0 / s->zoom_props->width[i]; } } /** * @brief Report Zoom Region Properties * * This function prints out a table containing the properties of the * zoom region, if it is enabled. The table includes information such as * dimensions, center, CDIM, background CDIM, buffer CDIM, region buffer * ratio, zoom boost factor, minimum zoom cell width, background cell width, * buffer width, and the number of wanderers. * * @param s The space */ void zoom_report_cell_properties(const struct space *s) { struct zoom_region_properties *zoom_props = s->zoom_props; /* Compute the number of background cells along each side of the void * region. */ int void_cdim[3]; for (int i = 0; i < 3; i++) { void_cdim[i] = (int)ceil(s->zoom_props->void_dim[i] * s->iwidth[i]); } /* Cdims */ message("%28s = [%d, %d, %d]", "Background cdim", s->cdim[0], s->cdim[1], s->cdim[2]); if (zoom_props->with_buffer_cells) message("%28s = [%d, %d, %d]", "Buffer cdim", zoom_props->buffer_cdim[0], zoom_props->buffer_cdim[1], zoom_props->buffer_cdim[2]); message("%28s = [%d, %d, %d]", "Zoom cdim", zoom_props->cdim[0], zoom_props->cdim[1], zoom_props->cdim[2]); message("%28s = [%d, %d, %d]", "Void cdim", void_cdim[0], void_cdim[1], void_cdim[2]); /* Dimensions */ message("%28s = [%f, %f, %f]", "Background Dimensions", s->dim[0], s->dim[1], s->dim[2]); if (zoom_props->with_buffer_cells) message("%28s = [%f, %f, %f]", "Buffer Region Dimensions", zoom_props->buffer_dim[0], zoom_props->buffer_dim[1], zoom_props->buffer_dim[2]); message("%28s = [%f, %f, %f]", "Zoom Region Dimensions", zoom_props->dim[0], zoom_props->dim[1], zoom_props->dim[2]); message("%28s = [%f, %f, %f]", "Void Region Dimensions", s->zoom_props->void_dim[0], s->zoom_props->void_dim[1], s->zoom_props->void_dim[2]); /* Cell Widths */ message("%28s = [%f, %f, %f]", "Background Cell Width", s->width[0], s->width[1], s->width[2]); if (zoom_props->with_buffer_cells) message("%28s = [%f, %f, %f]", "Buffer Cell Width", zoom_props->buffer_width[0], zoom_props->buffer_width[1], zoom_props->buffer_width[2]); message("%28s = [%f, %f, %f]", "Zoom Cell Width", zoom_props->width[0], zoom_props->width[1], zoom_props->width[2]); /* Number of Cells */ message("%28s = %d", "Number of Background Cells", zoom_props->nr_bkg_cells); if (zoom_props->with_buffer_cells) message("%28s = %d", "Number of Buffer Cells", zoom_props->nr_buffer_cells); message("%28s = %d", "Number of Zoom Cells", zoom_props->nr_zoom_cells); /* Bounds */ if (zoom_props->with_buffer_cells) message( "%28s = [%f-%f, %f-%f, %f-%f]", "Buffer Bounds", zoom_props->buffer_lower_bounds[0], zoom_props->buffer_upper_bounds[0], zoom_props->buffer_lower_bounds[1], zoom_props->buffer_upper_bounds[1], zoom_props->buffer_lower_bounds[2], zoom_props->buffer_upper_bounds[2]); message( "%28s = [%f-%f, %f-%f, %f-%f]", "Zoom Region Bounds", zoom_props->region_lower_bounds[0], zoom_props->region_upper_bounds[0], zoom_props->region_lower_bounds[1], zoom_props->region_upper_bounds[1], zoom_props->region_lower_bounds[2], zoom_props->region_upper_bounds[2]); /* Depths */ if (zoom_props->with_buffer_cells) message("%28s = %d", "Buffer Top Level Depth", zoom_props->buffer_cell_depth); message("%28s = %d", "Zoom Top Level Depth", zoom_props->zoom_cell_depth); message("%28s = %d", "Neighbour Max Tree Depth", zoom_props->neighbour_max_tree_depth); /* Assorted extra zoom properties */ message("%28s = %f", "Zoom Region Pad Factor", zoom_props->region_pad_factor); message("%28s = [%f, %f, %f]", "Zoom Region Shift", zoom_props->zoom_shift[0], zoom_props->zoom_shift[1], zoom_props->zoom_shift[2]); message("%28s = [%f, %f, %f]", "Zoom Region Center", zoom_props->region_lower_bounds[0] + (zoom_props->dim[0] / 2.0), zoom_props->region_lower_bounds[1] + (zoom_props->dim[1] / 2.0), zoom_props->region_lower_bounds[2] + (zoom_props->dim[2] / 2.0)); } /** * @brief Parse and set the zoom region properties. * * This function allocates the zoom region properties struct and populates it. * * If we're not running a zoom this function will do nothing. * * @param params Swift parameter structure. * @param s The space * @param verbose Are we talking? */ void zoom_props_init(struct swift_params *params, struct space *s, const int verbose) { /* If not, we're done here */ if (!s->with_zoom_region) { return; } /* Zoom region properties are stored in a structure. */ s->zoom_props = (struct zoom_region_properties *)malloc( sizeof(struct zoom_region_properties)); if (s->zoom_props == NULL) error("Error allocating memory for the zoom parameters."); bzero(s->zoom_props, sizeof(struct zoom_region_properties)); /* Calculate the gravity mesh distance, we need this for buffer cells and * neighbour cell labbeling later on. */ /* NOTE: when this is first called we don't have the gravity properties (and * the engine isn't attached to the space) yet so we need to read directly * from the params. */ /* Get the mesh size */ int mesh_size = parser_get_param_int(params, "Gravity:mesh_side_length"); /* Calculate the maximum distance at which we have a gravity task based * on the . */ float a_smooth = parser_get_opt_param_float(params, "Gravity:a_smooth", 1.25); float r_cut_max_ratio = parser_get_opt_param_float(params, "Gravity:r_cut_max", 4.5); float r_s = a_smooth * s->dim[0] / mesh_size; s->zoom_props->neighbour_distance = r_s * r_cut_max_ratio; /* Parse the parameter file and populate the properties struct. */ zoom_parse_params(params, s->zoom_props); } /** * @brief Initialise the zoom region geometry. * * This will compute the cell grid properties ready for cell * cosntruction when zoom_construct_tl_cells. * * @param s The space. * @param verbose Are we talking? */ void zoom_region_init(struct space *s, const int verbose) { /* Nothing to do if we are restarting, just report geometry and move on. */ if (s->e != NULL && s->e->restarting) { if (verbose) zoom_report_cell_properties(s); return; } /* Update the neighbour distance in case the gravity props have changed. */ if (s->e != NULL) { s->zoom_props->neighbour_distance = s->e->gravity_properties->r_s * s->e->gravity_properties->r_cut_max_ratio; } /* Compute the extent of the zoom region. * NOTE: this calculates the shift necessary to move the zoom region to * the centre of the box and stores it in s->zoom_props */ double ini_dim = zoom_get_region_dim_and_shift(s); /* Apply the shift to the particles. */ for (size_t k = 0; k < s->nr_parts; k++) { s->parts[k].x[0] += s->zoom_props->zoom_shift[0]; s->parts[k].x[1] += s->zoom_props->zoom_shift[1]; s->parts[k].x[2] += s->zoom_props->zoom_shift[2]; } for (size_t k = 0; k < s->nr_gparts; k++) { s->gparts[k].x[0] += s->zoom_props->zoom_shift[0]; s->gparts[k].x[1] += s->zoom_props->zoom_shift[1]; s->gparts[k].x[2] += s->zoom_props->zoom_shift[2]; } for (size_t k = 0; k < s->nr_sparts; k++) { s->sparts[k].x[0] += s->zoom_props->zoom_shift[0]; s->sparts[k].x[1] += s->zoom_props->zoom_shift[1]; s->sparts[k].x[2] += s->zoom_props->zoom_shift[2]; } for (size_t k = 0; k < s->nr_bparts; k++) { s->bparts[k].x[0] += s->zoom_props->zoom_shift[0]; s->bparts[k].x[1] += s->zoom_props->zoom_shift[1]; s->bparts[k].x[2] += s->zoom_props->zoom_shift[2]; } for (size_t k = 0; k < s->nr_sinks; k++) { s->sinks[k].x[0] += s->zoom_props->zoom_shift[0]; s->sinks[k].x[1] += s->zoom_props->zoom_shift[1]; s->sinks[k].x[2] += s->zoom_props->zoom_shift[2]; } /* Include the requested padding around the high resolution particles. */ double max_dim = ini_dim * s->zoom_props->region_pad_factor; /* Define the background grid (we'll treat this as gospel). */ for (int i = 0; i < 3; i++) { s->cdim[i] = s->zoom_props->bkg_cdim[i]; s->width[i] = s->dim[i] / s->cdim[i]; s->iwidth[i] = 1.0 / s->width[i]; } /* Compute the void region bounds and number of zoom regions that tesselate * it. */ int nr_zoom_regions = zoom_get_void_geometry(s, max_dim); /* Check the user gave a sensible background cdim, if the number of zoom * regions is too high we will have to set up a unworkable number of buffer * cells. */ if (nr_zoom_regions >= 64) { error( "Background cell size is too large relative to the zoom region! " "Increase ZoomRegion:bkg_top_level_cells (would have had %d zoom " "regions in the void region).", nr_zoom_regions); } /* If its alot but not silly just warn the user. */ if (nr_zoom_regions >= 16) { warning( "Background cell size is large relative to the zoom region! " "(would have had %d zoom regions in the void region). ", nr_zoom_regions); } if (verbose) { message("Initial geometry gives %d zoom regions in the void region.", nr_zoom_regions); } /* Construct the zoom region geometry. */ /* NOTE: here we entirely avoid any buffer cell considerations since they * provide no performance benefit and are for now vestigual. */ zoom_get_geometry_no_buffer_cells(s); /* Store what the true boost factor ended up being */ double input_pad_factor = s->zoom_props->region_pad_factor; s->zoom_props->region_pad_factor = s->zoom_props->dim[0] / ini_dim; /* Ensure we haven't got a zoom region smaller than the high resolution * particle distribution. */ if (s->zoom_props->dim[0] < ini_dim) { error( "Found a zoom region smaller than the high resolution particle " "distribution! Adjust the cell structure " "(ZoomRegion:bkg_top_level_cells, ZoomRegion:zoom_top_level_cells)"); } /* Let's be safe and warn if we have drastically changed the size of the * requested padding region. */ if ((s->zoom_props->region_pad_factor / input_pad_factor) >= 2) warning( "The pad region has to be %d times larger than requested. " "Either increase ZoomRegion:region_pad_factor, increase the " "number of background cells, or increase the depths of the zoom cells.", (int)(s->zoom_props->region_pad_factor / input_pad_factor)); /* If we didn't get an explicit neighbour cell depth we'll match the zoom * depth. */ s->zoom_props->neighbour_max_tree_depth = (s->zoom_props->neighbour_max_tree_depth < 0) ? s->zoom_props->zoom_cell_depth : s->zoom_props->neighbour_max_tree_depth; /* The neighbour depth must be less the zoom depth or higher if given. */ if (s->zoom_props->neighbour_max_tree_depth < s->zoom_props->zoom_cell_depth) { error( "Zoom neighbour cell depth (%d) must be greater than or equal to the " "zoom cell depth (%d).", s->zoom_props->neighbour_max_tree_depth, s->zoom_props->zoom_cell_depth); } /* Set the minimum allowed zoom cell width. */ const double zoom_dmax = max3(s->zoom_props->dim[0], s->zoom_props->dim[1], s->zoom_props->dim[2]); s->zoom_props->cell_min = 0.99 * zoom_dmax / s->zoom_props->cdim[0]; /* Set the minimum background cell size. */ const double dmax = max3(s->dim[0], s->dim[1], s->dim[2]); s->cell_min = 0.99 * dmax / s->cdim[0]; /* Store cell numbers and offsets. */ s->zoom_props->bkg_cell_offset = s->zoom_props->cdim[0] * s->zoom_props->cdim[1] * s->zoom_props->cdim[2]; s->zoom_props->nr_zoom_cells = s->zoom_props->bkg_cell_offset; s->zoom_props->nr_bkg_cells = s->cdim[0] * s->cdim[1] * s->cdim[2]; s->zoom_props->buffer_cell_offset = s->zoom_props->bkg_cell_offset + s->zoom_props->nr_bkg_cells; s->zoom_props->nr_buffer_cells = s->zoom_props->buffer_cdim[0] * s->zoom_props->buffer_cdim[1] * s->zoom_props->buffer_cdim[2]; /* Report what we have done */ if (verbose) { zoom_report_cell_properties(s); } }