Commit de07fd91 authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Merge branch 'minimal_smoothing_length' into 'master'

Implement an optional minimal smoothing length based on the softening of the particle

See merge request !732
parents 97f3111f 9d95c26c
......@@ -30,12 +30,16 @@
#include "dimension.h"
#include "equation_of_state.h"
#include "error.h"
#include "gravity_properties.h"
#include "hydro.h"
#include "kernel_hydro.h"
#include "parser.h"
#include "units.h"
#define hydro_props_default_max_iterations 30
#define hydro_props_default_volume_change 1.4f
#define hydro_props_default_h_max FLT_MAX
#define hydro_props_default_h_min_ratio 0.f
#define hydro_props_default_h_tolerance 1e-4
#define hydro_props_default_init_temp 0.f
#define hydro_props_default_min_temp 0.f
......@@ -104,6 +108,13 @@ void hydro_props_init(struct hydro_props *p,
p->h_max = parser_get_opt_param_float(params, "SPH:h_max",
hydro_props_default_h_max);
/* Minimal smoothing length ratio to softening */
p->h_min_ratio = parser_get_opt_param_float(params, "SPH:h_min_ratio",
hydro_props_default_h_min_ratio);
/* Temporarily set the minimal softening to 0. */
p->h_min = 0.f;
/* 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);
......@@ -317,6 +328,7 @@ void hydro_props_print_snapshot(hid_t h_grpsph, const struct hydro_props *p) {
* @param p the struct
*/
void hydro_props_init_no_hydro(struct hydro_props *p) {
p->eta_neighbours = 1.2348;
p->h_tolerance = hydro_props_default_h_tolerance;
p->target_neighbours = pow_dimension(p->eta_neighbours) * kernel_norm;
......@@ -325,6 +337,8 @@ void hydro_props_init_no_hydro(struct hydro_props *p) {
(pow_dimension(delta_eta) - pow_dimension(p->eta_neighbours)) *
kernel_norm;
p->h_max = hydro_props_default_h_max;
p->h_min = 0.f;
p->h_min_ratio = hydro_props_default_h_min_ratio;
p->max_smoothing_iterations = hydro_props_default_max_iterations;
p->CFL_condition = 0.1;
p->log_max_h_change = logf(powf(1.4, hydro_dimension_inv));
......@@ -352,6 +366,20 @@ void hydro_props_init_no_hydro(struct hydro_props *p) {
p->diffusion.alpha_min = hydro_props_default_diffusion_alpha_min;
}
/**
* @brief Update the global properties of the hydro scheme for that time-step.
*
* @param p The properties to update.
* @param gp The properties of the gravity scheme.
* @param cosmo The cosmological model.
*/
void hydro_props_update(struct hydro_props *p, const struct gravity_props *gp,
const struct cosmology *cosmo) {
/* Update the minimal allowed smoothing length */
p->h_min = p->h_min_ratio * gp->epsilon_cur;
}
/**
* @brief Write a hydro_props struct to the given FILE as a stream of bytes.
*
......
......@@ -32,10 +32,14 @@
#endif
/* Local includes. */
#include "parser.h"
#include "physical_constants.h"
#include "restart.h"
#include "units.h"
/* Forward declarations */
struct cosmology;
struct swift_params;
struct gravity_props;
struct phys_const;
struct unit_system;
/**
* @brief Contains all the constants and parameters of the hydro scheme
......@@ -57,6 +61,12 @@ struct hydro_props {
/*! Maximal smoothing length */
float h_max;
/*! Minimal smoothing length expressed as ratio to softening length */
float h_min_ratio;
/*! Minimal smoothing length */
float h_min;
/*! Maximal number of iterations to converge h */
int max_smoothing_iterations;
......@@ -131,6 +141,9 @@ void hydro_props_init(struct hydro_props *p,
const struct unit_system *us,
struct swift_params *params);
void hydro_props_update(struct hydro_props *p, const struct gravity_props *gp,
const struct cosmology *cosmo);
#if defined(HAVE_HDF5)
void hydro_props_print_snapshot(hid_t h_grpsph, const struct hydro_props *p);
#endif
......
......@@ -1230,6 +1230,7 @@ void runner_do_ghost(struct runner *r, struct cell *c, int timer) {
const struct cosmology *cosmo = e->cosmology;
const struct chemistry_global_data *chemistry = e->chemistry;
const float hydro_h_max = e->hydro_properties->h_max;
const float hydro_h_min = e->hydro_properties->h_min;
const float eps = e->hydro_properties->h_tolerance;
const float hydro_eta_dim =
pow_dimension(e->hydro_properties->eta_neighbours);
......@@ -1294,6 +1295,7 @@ void runner_do_ghost(struct runner *r, struct cell *c, int timer) {
const float h_old = p->h;
const float h_old_dim = pow_dimension(h_old);
const float h_old_dim_minus_one = pow_dimension_minus_one(h_old);
float h_new;
int has_no_neighbours = 0;
......@@ -1330,11 +1332,13 @@ void runner_do_ghost(struct runner *r, struct cell *c, int timer) {
#endif
/* Skip if h is already h_max and we don't have enough neighbours */
if ((p->h >= hydro_h_max) && (f < 0.f)) {
if (((p->h >= hydro_h_max) && (f < 0.f)) ||
((p->h <= hydro_h_min) && (f > 0.f))) {
/* We have a particle whose smoothing length is already set (wants
* to be larger but has already hit the maximum). So, just tidy up
* as if the smoothing length had converged correctly */
* to be larger but has already hit the maximum OR wants to be smaller
* but has already reached the minimum). So, just tidy up as if the
* smoothing length had converged correctly */
#ifdef EXTRA_HYDRO_LOOP
......@@ -1436,8 +1440,8 @@ void runner_do_ghost(struct runner *r, struct cell *c, int timer) {
p->h = h_new;
}
/* If below the absolute maximum, try again */
if (p->h < hydro_h_max) {
/* If within the allowed range, try again */
if (p->h < hydro_h_max && p->h > hydro_h_min) {
/* Flag for another round of fun */
pid[redo] = pid[i];
......@@ -1453,7 +1457,18 @@ void runner_do_ghost(struct runner *r, struct cell *c, int timer) {
/* Off we go ! */
continue;
} else {
} else if (p->h <= hydro_h_min) {
/* Ok, this particle is a lost cause... */
p->h = hydro_h_min;
/* Do some damage control if no neighbours at all were found */
if (has_no_neighbours) {
hydro_part_has_no_neighbours(p, xp, cosmo);
chemistry_part_has_no_neighbours(p, xp, chemistry, cosmo);
}
} else if (p->h >= hydro_h_max) {
/* Ok, this particle is a lost cause... */
p->h = hydro_h_max;
......
......@@ -3351,6 +3351,10 @@ void space_first_init_parts_mapper(void *restrict map_data, int count,
const struct hydro_props *hydro_props = s->e->hydro_properties;
const float u_init = hydro_props->initial_internal_energy;
const float hydro_h_min_ratio = e->hydro_properties->h_min_ratio;
const struct gravity_props *grav_props = s->e->gravity_properties;
const int with_gravity = e->policy & engine_policy_self_gravity;
const struct chemistry_global_data *chemistry = e->chemistry;
const struct cooling_function_data *cool_func = e->cooling_func;
......@@ -3360,6 +3364,12 @@ void space_first_init_parts_mapper(void *restrict map_data, int count,
if (p[k].h <= 0.)
error("Invalid value of smoothing length for part %lld h=%e", p[k].id,
p[k].h);
if (with_gravity) {
const struct gpart *gp = p[k].gpart;
const float softening = gravity_get_softening(gp, grav_props);
p->h = max(p->h, softening * hydro_h_min_ratio);
}
}
/* Convert velocities to internal units */
......
Markdown is supported
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