Commit 37baeabb authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Prevent the user from running with a mesh that is too small given r_cut_max.

parent 4e8b4efd
...@@ -58,12 +58,17 @@ void gravity_props_init(struct gravity_props *p, struct swift_params *params, ...@@ -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( p->r_cut_min_ratio = parser_get_opt_param_float(
params, "Gravity:r_cut_min", gravity_props_default_r_cut_min); params, "Gravity:r_cut_min", gravity_props_default_r_cut_min);
/* Some basic checks */
if (p->mesh_size % 2 != 0) if (p->mesh_size % 2 != 0)
error("The mesh side-length must be an even number."); error("The mesh side-length must be an even number.");
if (p->a_smooth <= 0.) if (p->a_smooth <= 0.)
error("The mesh smoothing scale 'a_smooth' must be > 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 */ /* Time integration */
p->eta = parser_get_param_float(params, "Gravity:eta"); p->eta = parser_get_param_float(params, "Gravity:eta");
......
...@@ -574,7 +574,7 @@ void pm_mesh_init(struct pm_mesh* mesh, const struct gravity_props* props, ...@@ -574,7 +574,7 @@ void pm_mesh_init(struct pm_mesh* mesh, const struct gravity_props* props,
#ifdef HAVE_FFTW #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"); error("Doing mesh-gravity on a non-cubic domain");
const int N = props->mesh_size; const int N = props->mesh_size;
...@@ -591,6 +591,9 @@ void pm_mesh_init(struct pm_mesh* mesh, const struct gravity_props* props, ...@@ -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_max = mesh->r_s * props->r_cut_max_ratio;
mesh->r_cut_min = mesh->r_s * props->r_cut_min_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 */ /* Allocate the memory for the combined density and potential array */
mesh->potential = (double*)fftw_malloc(sizeof(double) * N * N * N); mesh->potential = (double*)fftw_malloc(sizeof(double) * N * N * N);
if (mesh->potential == NULL) if (mesh->potential == NULL)
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment