From 391c962b84acfe1a3c08f83fb6cd916fc7cbf620 Mon Sep 17 00:00:00 2001
From: Matthieu Schaller <matthieu.schaller@durham.ac.uk>
Date: Fri, 3 Mar 2017 16:00:57 +0000
Subject: [PATCH] Added a maximal smoothing length that can be used to prevent
particles from becoming too big in vacuum cases near the box edges. Will
break the accuracy of the hydrodynamic scheme but is necessary in some
situations.
---
examples/parameter_example.yml | 2 +-
src/cell.c | 18 +++++++++++-------
src/hydro_properties.c | 11 ++++++++++-
src/hydro_properties.h | 5 ++++-
src/runner.c | 26 +++++++++++++++++---------
5 files changed, 43 insertions(+), 19 deletions(-)
diff --git a/examples/parameter_example.yml b/examples/parameter_example.yml
index bdcd880c83..7c69d19bc2 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 517386a460..ab5b99c030 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 8db23a3de0..46785b4b2d 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 859ae10f18..716c4c060c 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 2e581e029f..dc83185593 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 */
--
GitLab