Commit cf3d11cf authored by Matthieu Schaller's avatar Matthieu Schaller

Updated the way we converge towards the smoothing length. We now use a...

Updated the way we converge towards the smoothing length. We now use a Newton-Raphson scheme. This applies to the Gadget-2 scheme.
parent fda9dfdc
......@@ -43,11 +43,11 @@ Statistics:
# Parameters for the hydrodynamics scheme
SPH:
resolution_eta: 1.2348 # Target smoothing length in units of the mean inter-particle separation (1.2348 == 48Ngbs with the cubic spline kernel).
delta_neighbours: 0.1 # The tolerance for the targetted number of neighbours.
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: 1.4 # (Optional) Maximal allowed change of kernel volume over one time-step
h_tolerance: 1e-4 # (Optional) Relative accuracy of the Netwon-Raphson scheme for the smoothing lengths.
h_max: 10. # (Optional) Maximal allowed smoothing length in internal units. Defaults to FLT_MAX if unspecified.
max_volume_change: 1.4 # (Optional) Maximal allowed change of kernel volume over one time-step.
max_ghost_iterations: 30 # (Optional) Maximal number of iterations allowed to converge towards the smoothing length.
# Parameters for the self-gravity scheme
Gravity:
......
......@@ -1718,7 +1718,7 @@ int cell_unskip_tasks(struct cell *c, struct scheduler *s) {
else
error("Drift task missing !");
if (t->type == task_type_pair) scheduler_activate(s, cj->drift_part);
if (cell_is_active(cj, e)) {
for (l = cj->send_rho; l != NULL && l->t->cj->nodeID != ci->nodeID;
......@@ -1870,11 +1870,10 @@ void cell_drift_part(struct cell *c, const struct engine *e) {
float dx_max = 0.f, dx2_max = 0.f;
float dx_max_sort = 0.0f, dx2_max_sort = 0.f;
float cell_h_max = 0.f;
#ifdef SWIFT_DEBUG_CHECKS
/* Check that we only drift local cells. */
if (c->nodeID != engine_rank)
error("Drifting a foreign cell is nope.");
if (c->nodeID != engine_rank) error("Drifting a foreign cell is nope.");
/* Check that we are actually going to move forward. */
if (ti_current < ti_old_part) error("Attempt to drift to the past");
......
......@@ -118,6 +118,34 @@ __attribute__((always_inline)) INLINE static float pow_dimension_plus_one(
#endif
}
/**
* @brief Returns the argument to the power given by the dimension minus one
*
* Computes \f$x^{d-1}\f$.
*/
__attribute__((always_inline)) INLINE static float pow_dimension_minus_one(
float x) {
#if defined(HYDRO_DIMENSION_3D)
return x * x;
#elif defined(HYDRO_DIMENSION_2D)
return x;
#elif defined(HYDRO_DIMENSION_1D)
return 1.f;
#else
error("The dimension is not defined !");
return 0.f;
#endif
}
/**
* @brief Inverts the given dimension by dimension matrix (in place)
*
......
......@@ -4266,7 +4266,7 @@ void engine_init(struct engine *e, struct space *s,
"Version: %s \n# "
"Number of threads: %d\n# Number of MPI ranks: %d\n# Hydrodynamic "
"scheme: %s\n# Hydrodynamic kernel: %s\n# No. of neighbours: %.2f "
"+/- %.2f\n# Eta: %f\n",
"+/- %.4f\n# Eta: %f\n",
hostname(), git_branch(), git_revision(), compiler_name(),
compiler_version(), e->nr_threads, e->nr_nodes, SPH_IMPLEMENTATION,
kernel_name, e->hydro_properties->target_neighbours,
......
......@@ -206,12 +206,13 @@ __attribute__((always_inline)) INLINE static void hydro_end_density(
p->rho += p->mass * kernel_root;
p->density.rho_dh -= hydro_dimension * p->mass * kernel_root;
p->density.wcount += kernel_root;
p->density.wcount_dh -= hydro_dimension * kernel_root;
/* Finish the calculation by inserting the missing h-factors */
p->rho *= h_inv_dim;
p->density.rho_dh *= h_inv_dim_plus_one;
p->density.wcount *= kernel_norm;
p->density.wcount_dh *= h_inv * kernel_gamma * kernel_norm;
p->density.wcount *= h_inv_dim;
p->density.wcount_dh *= h_inv_dim_plus_one;
const float rho_inv = 1.f / p->rho;
......@@ -224,6 +225,31 @@ __attribute__((always_inline)) INLINE static void hydro_end_density(
p->density.div_v *= h_inv_dim_plus_one * rho_inv;
}
/**
* @brief Sets all particle fields to sensible values when the #part has 0 ngbs.
*
* @param p The particle to act upon
* @param xp The extended particle data to act upon
*/
__attribute__((always_inline)) INLINE static void hydro_part_has_no_neighbours(
struct part *restrict p, struct xpart *restrict xp) {
/* Some smoothing length multiples. */
const float h = p->h;
const float h_inv = 1.0f / h; /* 1/h */
const float h_inv_dim = pow_dimension(h_inv); /* 1/h^d */
/* Re-set problematic values */
p->rho = p->mass * kernel_root * h_inv_dim;
p->density.wcount = kernel_root * kernel_norm * h_inv_dim;
p->density.rho_dh = 0.f;
p->density.wcount_dh = 0.f;
p->density.div_v = 0.f;
p->density.rot_v[0] = 0.f;
p->density.rot_v[1] = 0.f;
p->density.rot_v[2] = 0.f;
}
/**
* @brief Prepare a particle for the force calculation.
*
......@@ -239,6 +265,9 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
const float fac_mu = 1.f; /* Will change with cosmological integration */
/* Inverse of the physical density */
const float rho_inv = 1.f / p->rho;
/* Compute the norm of the curl */
const float curl_v = sqrtf(p->density.rot_v[0] * p->density.rot_v[0] +
p->density.rot_v[1] * p->density.rot_v[1] +
......@@ -254,7 +283,6 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
const float soundspeed = gas_soundspeed_from_pressure(p->rho, pressure);
/* Divide the pressure by the density squared to get the SPH term */
const float rho_inv = 1.f / p->rho;
const float P_over_rho2 = pressure * rho_inv * rho_inv;
/* Compute the Balsara switch */
......@@ -262,11 +290,11 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
abs_div_v / (abs_div_v + curl_v + 0.0001f * soundspeed / fac_mu / p->h);
/* Compute the "grad h" term */
const float grad_h_term =
const float omega_inv =
1.f / (1.f + hydro_dimension_inv * p->h * p->density.rho_dh * rho_inv);
/* Update variables. */
p->force.f = grad_h_term;
p->force.f = omega_inv;
p->force.P_over_rho2 = P_over_rho2;
p->force.soundspeed = soundspeed;
p->force.balsara = balsara;
......
This diff is collapsed.
......@@ -35,14 +35,24 @@
#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_tolerance 1e-4
void hydro_props_init(struct hydro_props *p,
const struct swift_params *params) {
/* Kernel properties */
p->eta_neighbours = parser_get_param_float(params, "SPH:resolution_eta");
/* Tolerance for the smoothing length Newton-Raphson scheme */
p->h_tolerance = parser_get_opt_param_float(params, "SPH:h_tolerance",
hydro_props_default_h_tolerance);
/* Get derived properties */
p->target_neighbours = pow_dimension(p->eta_neighbours) * kernel_norm;
p->delta_neighbours = parser_get_param_float(params, "SPH:delta_neighbours");
const float delta_eta = p->eta_neighbours * (1.f + p->h_tolerance);
p->delta_neighbours =
(pow_dimension(delta_eta) - pow_dimension(p->eta_neighbours)) *
kernel_norm;
#ifdef SHADOWFAX_SPH
/* change the meaning of target_neighbours and delta_neighbours */
......@@ -81,9 +91,11 @@ void hydro_props_print(const struct hydro_props *p) {
message("Hydrodynamic scheme: %s in %dD.", SPH_IMPLEMENTATION,
(int)hydro_dimension);
message("Hydrodynamic kernel: %s with %.2f +/- %.2f neighbours (eta=%f).",
kernel_name, p->target_neighbours, p->delta_neighbours,
p->eta_neighbours);
message("Hydrodynamic kernel: %s with eta=%f (%.2f neighbours).", kernel_name,
p->eta_neighbours, p->target_neighbours);
message("Hydrodynamic tolerance in h: %.5f (+/- %.4f neighbours).",
p->h_tolerance, p->delta_neighbours);
message("Hydrodynamic integration: CFL parameter: %.4f.", p->CFL_condition);
......@@ -110,6 +122,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, "Smoothing length tolerance", p->h_tolerance);
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))",
......
......@@ -16,10 +16,14 @@
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*
******************************************************************************/
#ifndef SWIFT_HYDRO_PROPERTIES
#define SWIFT_HYDRO_PROPERTIES
/**
* @file hydro_properties.h
* @brief Contains all the constants and parameters of the hydro scheme
*/
/* Config parameters. */
#include "../config.h"
......@@ -35,19 +39,28 @@
*/
struct hydro_props {
/* Kernel properties */
/*! Resolution parameter */
float eta_neighbours;
/*! Target weightd number of neighbours (for info only)*/
float target_neighbours;
/*! Smoothing length tolerance */
float h_tolerance;
/*! Tolerance on neighbour number (for info only)*/
float delta_neighbours;
/* Maximal smoothing length */
/*! Maximal smoothing length */
float h_max;
/* Number of iterations to converge h */
/*! Maximal number of iterations to converge h */
int max_smoothing_iterations;
/* Time integration properties */
/*! Time integration properties */
float CFL_condition;
/*! Maximal change of h over one time-step */
float log_max_h_change;
};
......
......@@ -622,11 +622,9 @@ void runner_do_ghost(struct runner *r, struct cell *c, int timer) {
const struct engine *e = r->e;
const struct space *s = e->s;
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;
const float min_wcount =
target_wcount - e->hydro_properties->delta_neighbours;
const float eps = e->hydro_properties->h_tolerance;
const float hydro_eta_dim =
pow_dimension(e->hydro_properties->eta_neighbours);
const int max_smoothing_iter = e->hydro_properties->max_smoothing_iterations;
int redo = 0, count = 0;
......@@ -670,28 +668,47 @@ void runner_do_ghost(struct runner *r, struct cell *c, int timer) {
if (!part_is_active(p, e)) error("Ghost applied to inactive particle");
#endif
/* Finish the density calculation */
hydro_end_density(p);
/* Get some useful values */
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;
/* Did we get the right number of neighbours? */
if (p->density.wcount > max_wcount || p->density.wcount < min_wcount) {
if (p->density.wcount == 0.f) { /* No neighbours case */
float h_corr = 0.f;
/* Double h and try again */
h_new = 2.f * h_old;
} else {
/* If no derivative, double the smoothing length. */
if (p->density.wcount_dh == 0.0f) h_corr = p->h;
/* Finish the density calculation */
hydro_end_density(p);
/* Otherwise, compute the smoothing length update (Newton step). */
else {
h_corr = (target_wcount - p->density.wcount) / p->density.wcount_dh;
/* Compute one step of the Newton-Raphson scheme */
const float n_sum = p->density.wcount * h_old_dim;
const float n_target = hydro_eta_dim;
const float f = n_sum - n_target;
const float f_prime =
p->density.wcount_dh * h_old_dim +
hydro_dimension * p->density.wcount * h_old_dim_minus_one;
/* Truncate to the range [ -p->h/2 , p->h ]. */
h_corr = (h_corr < p->h) ? h_corr : p->h;
h_corr = (h_corr > -0.5f * p->h) ? h_corr : -0.5f * p->h;
}
h_new = h_old - f / f_prime;
#ifdef SWIFT_DEBUG_CHECKS
if ((f > 0.f && h_new > h_old) || (f < 0.f && h_new < h_old))
error(
"Smoothing length correction not going in the right direction");
#endif
/* Safety check: truncate to the range [ h_old/2 , 2h_old ]. */
h_new = min(h_new, 2.f * h_old);
h_new = max(h_new, 0.5f * h_old);
}
/* Check whether the particle has an inappropriate smoothing length */
if (fabsf(h_new - h_old) > eps * h_old) {
/* Ok, correct then */
p->h += h_corr;
p->h = h_new;
/* If below the absolute maximum, try again */
if (p->h < hydro_h_max) {
......@@ -709,6 +726,10 @@ void runner_do_ghost(struct runner *r, struct cell *c, int timer) {
/* Ok, this particle is a lost cause... */
p->h = hydro_h_max;
/* Do some damage control if no neighbours at all were found */
if (p->density.wcount == kernel_root * kernel_norm)
hydro_part_has_no_neighbours(p, xp);
}
}
......
......@@ -549,8 +549,7 @@ int main(int argc, char *argv[]) {
prog_const.const_newton_G = 1.f;
struct hydro_props hp;
hp.target_neighbours = pow_dimension(h) * kernel_norm;
hp.delta_neighbours = 4.;
hp.h_tolerance = 1e0;
hp.h_max = FLT_MAX;
hp.max_smoothing_iterations = 1;
hp.CFL_condition = 0.1;
......
# ID pos_x pos_y pos_z v_x v_y v_z rho rho_dh wcount wcount_dh div_v curl_vx curl_vy curl_vz
0 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 2e-6 4e-5 2e-4 2e-3 1e-5 6e-6 6e-6 6e-6
0 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 2e-6 4e-5 3e-4 1e-2 1e-5 6e-6 6e-6 6e-6
0 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1.2e-4 1e-4 1e-4 2e-4 1e-4 1e-4 1e-4
0 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6
# ID pos_x pos_y pos_z v_x v_y v_z rho rho_dh wcount wcount_dh div_v curl_vx curl_vy curl_vz
0 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 2e-6 1e-4 2e-4 3e-3 1e-5 3e-6 3e-6 7e-6
0 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 2e-3 1e-5 1e-4 6e-5 2e-3 2e-3 2e-3
0 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 4e-4 1e-6 1e-6 1e-6 2e-6 2e-6 2e-6
0 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 2e-6 1e-4 2e-4 1e-2 1e-5 3e-6 3e-6 7e-6
0 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 2e-3 1e-5 1e-3 6e-5 2e-3 2e-3 2e-3
0 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 4e-4 1e-6 1e0 1e-6 2e-6 2e-6 2e-6
# ID pos_x pos_y pos_z v_x v_y v_z rho rho_dh wcount wcount_dh div_v curl_vx curl_vy curl_vz
0 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 2e-6 1e-4 2e-4 4e-3 1e-5 3e-6 3e-6 7e-6
0 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1.4e-2 1e-5 1e-4 2.5e-4 3e-3 3e-3 3e-3
0 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 4e-6 4e-6 4e-6
0 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 2e-6 1e-4 4e-4 1.2e-2 1e-5 3e-6 3e-6 7e-6
0 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1.4e-2 1e-5 1e-3 2.5e-4 3e-3 3e-3 3e-3
0 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1e0 1e-6 4e-6 4e-6 4e-6
# ID pos_x pos_y pos_z v_x v_y v_z rho rho_dh wcount wcount_dh div_v curl_vx curl_vy curl_vz
0 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 4e-6 4e-5 3e-4 3e-3 2e-4 2e-4 2e-4 2e-4
0 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 2e-6 2e-4 1e-4 1e-4 6e-4 2e-3 2e-3 2e-3
0 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 4e-6 4e-5 1e-3 1e-2 2e-4 2e-4 2e-4 2e-4
0 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 2e-6 2e-4 1e-4 2e-4 6e-4 2e-3 2e-3 2e-3
0 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1e-4 1e-6 1e-4 5e-4 2e-4 2e-4 2e-4
# ID pos_x pos_y pos_z v_x v_y v_z rho rho_dh wcount wcount_dh div_v curl_vx curl_vy curl_vz
0 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 3e-6 4e-5 2e-4 3e-3 2e-4 1e-4 1e-4 1e-4
0 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 2e-6 6e-3 1e-4 1e-4 1e-2 6e-3 6e-3 6e-3
0 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 2e-3 1e-6 1e-4 5e-4 3e-3 3e-3 3e-3
0 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 3e-6 4e-5 1e-3 1e-2 2e-4 1e-4 1e-4 1e-4
0 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 2e-6 6e-3 1e-4 3e-3 1e-2 6e-3 6e-3 6e-3
0 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 2e-3 1e-6 1e0 5e-4 3e-3 3e-3 3e-3
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