diff --git a/src/gravity_properties.c b/src/gravity_properties.c index fc1ce1d62e02c32d44667d602448fc4eb3a65344..40cd7cecb10eab7ebbe7b59d34fb15dccac8b620 100644 --- a/src/gravity_properties.c +++ b/src/gravity_properties.c @@ -58,12 +58,17 @@ void gravity_props_init(struct gravity_props *p, struct swift_params *params, p->r_cut_min_ratio = parser_get_opt_param_float( params, "Gravity:r_cut_min", gravity_props_default_r_cut_min); + /* Some basic checks */ if (p->mesh_size % 2 != 0) error("The mesh side-length must be an even number."); if (p->a_smooth <= 0.) error("The mesh smoothing scale 'a_smooth' must be > 0."); + if (2. * p->a_smooth * p->r_cut_max_ratio > p->mesh_size) + error("Mesh too small given r_cut_max. Should be at least %d cells wide.", + (int)(2. * p->a_smooth * p->r_cut_max_ratio) + 1); + /* Time integration */ p->eta = parser_get_param_float(params, "Gravity:eta"); diff --git a/src/mesh_gravity.c b/src/mesh_gravity.c index ff0c3eb5ecc6a50b7ac9121f3c0179bcbf9bfcf5..68b08c26d6d4cc37400637a05f6166607812f3d7 100644 --- a/src/mesh_gravity.c +++ b/src/mesh_gravity.c @@ -574,7 +574,7 @@ void pm_mesh_init(struct pm_mesh* mesh, const struct gravity_props* props, #ifdef HAVE_FFTW - if (dim[0] != dim[1] || dim[0] != dim[1]) + if (dim[0] != dim[1] || dim[0] != dim[2]) error("Doing mesh-gravity on a non-cubic domain"); const int N = props->mesh_size; @@ -591,6 +591,9 @@ void pm_mesh_init(struct pm_mesh* mesh, const struct gravity_props* props, mesh->r_cut_max = mesh->r_s * props->r_cut_max_ratio; mesh->r_cut_min = mesh->r_s * props->r_cut_min_ratio; + if (2. * mesh->r_cut_max > box_size) + error("Mesh too small or r_cut_max too big for this box size"); + /* Allocate the memory for the combined density and potential array */ mesh->potential = (double*)fftw_malloc(sizeof(double) * N * N * N); if (mesh->potential == NULL)