Commit 391c962b authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Added a maximal smoothing length that can be used to prevent particles from...

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.
parent 2e7b35b1
......@@ -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:
......
......@@ -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 */
......
......@@ -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);
......
......@@ -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 */
......
......@@ -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 */
......
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