diff --git a/examples/parameter_example.yml b/examples/parameter_example.yml
index bdcd880c83eb02fdbed8d9d4da8642f3a90fc97b..7c69d19bc215472c90bb1b91f23f46702afc64f2 100644
--- a/examples/parameter_example.yml
+++ b/examples/parameter_example.yml
@@ -47,6 +47,7 @@ SPH:
CFL_condition: 0.1 # Courant-Friedrich-Levy condition for time integration.
max_ghost_iterations: 30 # (Optional) Maximal number of iterations allowed to converge towards the smoothing length.
max_volume_change: 2. # (Optional) Maximal allowed change of kernel volume over one time-step
+ h_max: 10. # (Optional) Maximal allowed smoothing length in internal units. Defaults to FLT_MAX if unspecified.
# Parameters for the self-gravity scheme
Gravity:
@@ -54,7 +55,6 @@ Gravity:
epsilon: 0.1 # Softening length (in internal units).
a_smooth: 1.25 # (Optional) Smoothing scale in top-level cell sizes to smooth the long-range forces over (this is the default value).
r_cut: 4.5 # (Optional) Cut-off in number of top-level cells beyond which no FMM forces are computed (this is the default value).
-
# Parameters related to the initial conditions
InitialConditions:
diff --git a/src/cell.c b/src/cell.c
index 517386a4608543ef0558e902107d931ccdeb0695..ab5b99c03024b1cc218c358b6d7051632496c665 100644
--- a/src/cell.c
+++ b/src/cell.c
@@ -1357,6 +1357,7 @@ void cell_set_super(struct cell *c, struct cell *super) {
*/
void cell_drift_particles(struct cell *c, const struct engine *e) {
+ const float hydro_h_max = e->hydro_properties->h_max;
const double timeBase = e->timeBase;
const integertime_t ti_old = c->ti_old;
const integertime_t ti_current = e->ti_current;
@@ -1367,7 +1368,7 @@ void cell_drift_particles(struct cell *c, const struct engine *e) {
/* Drift from the last time the cell was drifted to the current time */
const double dt = (ti_current - ti_old) * timeBase;
- float dx_max = 0.f, dx2_max = 0.f, h_max = 0.f;
+ float dx_max = 0.f, dx2_max = 0.f, cell_h_max = 0.f;
/* Check that we are actually going to move forward. */
if (ti_current < ti_old) error("Attempt to drift to the past");
@@ -1381,7 +1382,7 @@ void cell_drift_particles(struct cell *c, const struct engine *e) {
struct cell *cp = c->progeny[k];
cell_drift_particles(cp, e);
dx_max = max(dx_max, cp->dx_max);
- h_max = max(h_max, cp->h_max);
+ cell_h_max = max(cell_h_max, cp->h_max);
}
} else if (ti_current > ti_old) {
@@ -1400,7 +1401,7 @@ void cell_drift_particles(struct cell *c, const struct engine *e) {
const float dx2 = gp->x_diff[0] * gp->x_diff[0] +
gp->x_diff[1] * gp->x_diff[1] +
gp->x_diff[2] * gp->x_diff[2];
- dx2_max = (dx2_max > dx2) ? dx2_max : dx2;
+ dx2_max = max(dx2_max, dx2);
}
/* Loop over all the gas particles in the cell */
@@ -1414,14 +1415,17 @@ void cell_drift_particles(struct cell *c, const struct engine *e) {
/* Drift... */
drift_part(p, xp, dt, timeBase, ti_old, ti_current);
+ /* Limit h to within the allowed range */
+ p->h = min(p->h, hydro_h_max);
+
/* Compute (square of) motion since last cell construction */
const float dx2 = xp->x_diff[0] * xp->x_diff[0] +
xp->x_diff[1] * xp->x_diff[1] +
xp->x_diff[2] * xp->x_diff[2];
- dx2_max = (dx2_max > dx2) ? dx2_max : dx2;
+ dx2_max = max(dx2_max, dx2);
/* Maximal smoothing length */
- h_max = (h_max > p->h) ? h_max : p->h;
+ cell_h_max = max(cell_h_max, p->h);
}
/* Loop over all the star particles in the cell */
@@ -1442,12 +1446,12 @@ void cell_drift_particles(struct cell *c, const struct engine *e) {
} else {
- h_max = c->h_max;
+ cell_h_max = c->h_max;
dx_max = c->dx_max;
}
/* Store the values */
- c->h_max = h_max;
+ c->h_max = cell_h_max;
c->dx_max = dx_max;
/* Update the time of the last drift */
diff --git a/src/hydro_properties.c b/src/hydro_properties.c
index 8db23a3de0123400fe691b0d60f076929e1b0b48..46785b4b2d5b958f6db3bd9813d139575217d6fe 100644
--- a/src/hydro_properties.c
+++ b/src/hydro_properties.c
@@ -34,6 +34,7 @@
#define hydro_props_default_max_iterations 30
#define hydro_props_default_volume_change 2.0f
+#define hydro_props_default_h_max FLT_MAX
void hydro_props_init(struct hydro_props *p,
const struct swift_params *params) {
@@ -43,7 +44,11 @@ void hydro_props_init(struct hydro_props *p,
p->target_neighbours = pow_dimension(p->eta_neighbours) * kernel_norm;
p->delta_neighbours = parser_get_param_float(params, "SPH:delta_neighbours");
- /* Ghost stuff */
+ /* Maximal smoothing length */
+ p->h_max = parser_get_opt_param_float(params, "SPH:h_max",
+ hydro_props_default_h_max);
+
+ /* Number of iterations to converge h */
p->max_smoothing_iterations = parser_get_opt_param_int(
params, "SPH:max_ghost_iterations", hydro_props_default_max_iterations);
@@ -81,6 +86,9 @@ void hydro_props_print(const struct hydro_props *p) {
"(max|dlog(h)/dt|=%f).",
pow_dimension(expf(p->log_max_h_change)), p->log_max_h_change);
+ if (p->h_max != hydro_props_default_h_max)
+ message("Maximal smoothing length allowed: %.4f", p->h_max);
+
if (p->max_smoothing_iterations != hydro_props_default_max_iterations)
message("Maximal iterations in ghost task set to %d (default is %d)",
p->max_smoothing_iterations, hydro_props_default_max_iterations);
@@ -96,6 +104,7 @@ void hydro_props_print_snapshot(hid_t h_grpsph, const struct hydro_props *p) {
io_write_attribute_f(h_grpsph, "Kernel target N_ngb", p->target_neighbours);
io_write_attribute_f(h_grpsph, "Kernel delta N_ngb", p->delta_neighbours);
io_write_attribute_f(h_grpsph, "Kernel eta", p->eta_neighbours);
+ io_write_attribute_f(h_grpsph, "Maximal smoothing length", p->h_max);
io_write_attribute_f(h_grpsph, "CFL parameter", p->CFL_condition);
io_write_attribute_f(h_grpsph, "Volume log(max(delta h))",
p->log_max_h_change);
diff --git a/src/hydro_properties.h b/src/hydro_properties.h
index 859ae10f18a71ba920f44b8c59e63abd57cd7c92..716c4c060c21eb95d05f9d50e13d4681a958a6fd 100644
--- a/src/hydro_properties.h
+++ b/src/hydro_properties.h
@@ -40,7 +40,10 @@ struct hydro_props {
float target_neighbours;
float delta_neighbours;
- /* Kernel properties */
+ /* Maximal smoothing length */
+ float h_max;
+
+ /* Number of iterations to converge h */
int max_smoothing_iterations;
/* Time integration properties */
diff --git a/src/runner.c b/src/runner.c
index 2e581e029f4a5acbf6d9136c9f96925906064a5e..dc831855930f01dd745273e553e8f593c6238b92 100644
--- a/src/runner.c
+++ b/src/runner.c
@@ -496,8 +496,7 @@ void runner_do_init(struct runner *r, struct cell *c, int timer) {
if (!cell_is_active(c, e)) return;
/* Reset the gravity acceleration tensors */
- if (e->policy & engine_policy_self_gravity)
- multipole_init(c->multipole);
+ if (e->policy & engine_policy_self_gravity) multipole_init(c->multipole);
/* Recurse? */
if (c->split) {
@@ -596,6 +595,7 @@ void runner_do_ghost(struct runner *r, struct cell *c, int timer) {
struct part *restrict parts = c->parts;
struct xpart *restrict xparts = c->xparts;
const struct engine *e = r->e;
+ const float hydro_h_max = e->hydro_properties->h_max;
const float target_wcount = e->hydro_properties->target_neighbours;
const float max_wcount =
target_wcount + e->hydro_properties->delta_neighbours;
@@ -667,15 +667,23 @@ void runner_do_ghost(struct runner *r, struct cell *c, int timer) {
/* Ok, correct then */
p->h += h_corr;
- /* Flag for another round of fun */
- pid[redo] = pid[i];
- redo += 1;
+ /* If below the absolute maximum, try again */
+ if (p->h < hydro_h_max) {
- /* Re-initialise everything */
- hydro_init_part(p);
+ /* Flag for another round of fun */
+ pid[redo] = pid[i];
+ redo += 1;
- /* Off we go ! */
- continue;
+ /* Re-initialise everything */
+ hydro_init_part(p);
+
+ /* Off we go ! */
+ continue;
+ } else {
+
+ /* Ok, this particle is a lost cause... */
+ p->h = hydro_h_max;
+ }
}
/* We now have a particle whose smoothing length has converged */